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On the Stability of Two-Dimensional 
Smooth Transonic Flows 


YUNG-HUVAI KUO 
Cornell University 


SUMMARY 


Steady two-dimensional potential flows involving imbedded 
supercritical regions adjacent to a body have been calculated, but 
their existence has not yet been established experimentally. To 
further experimentation, the question of stability of such flows is 
examined. 

Following the standard procedure, a theoretically possible flow 
is disturbed at a certain instant, and the behavior of the subse- 
quent disturbance motion is investigated. By means of the 
transonic approximations, a simple solution valid in the neigh- 
borhood of the sonic line has been obtained. This study leads 
to the conclusion that (1) smooth decelerating transonic flows 
over bodies at negative curvature are unstable to compression 
pulses; (2) smooth accelerating transonic flows over any sur 
face are stable; and (3) smooth transonic flows with local super 


sonic region are unstable. 


(1) INTRODUCTION 


* that steady 


T HAS BEEN SHOWN THEORETICALLY! 
two-dimensional potential flows of gas past bodies 


may become locally supersonic at sufficiently high free- 
stream Mach Numbers. The supersonic regions are 
always located at the surface of a body, and their ex- 
tent increases as the Mach Number is raised. It is 
found, moreover, that the transition from subsonic to 
supersonic flow and vice versa takes place continually, 
and it appears that such smooth flows occur over a 
range of Mach Numbers which is terminated by the 
appearance of limiting lines. On the other hand, 
wind-tunnel observations show that such smooth trans- 
onic flows are generally replaced by flows involving 
Shock waves. The understanding and interpretation of 
this discrepancy between theory and experiment seems 

Presented at the Aerodynamics Session, Eighteenth Annual 
Meeting, I.A.S., New York, January 23-26, 1950. Revised 
and received September 15, 1950 P 

* This study was initiated by the work of A. Kantrowitz,!° 
who considered the stability of one-dimensional channel flow, and 
is being carried out with the financial assistance of the National 
The Committee has 


generously approved the publication of this preliminary report. 


Advisory Committee for Aeronautics. 


highly desirable and might lead to some device whereby 
the drag divergence and other undesirable phenomena 
associated with shock waves might be prevented or, 
at least, delayed to higher Mach Numbers. Hence, 
the establishment of criteria under which such flows are 
physically possible is not merely of academic interest 
but also of practical importance. 


Since, theoretically, the smooth transonic flow could 
exist for a much wider range of Mach Numbers than it 
actually does, the question of its stability arises. To 
answer this question, one must introduce disturbances 
at a certain moment and follow their future course. 
If the original steady motion is unstable, the disturb- 
ances will grow and eventually become shock waves, 
resulting in increase of drag. On the other hand, if it 
is stable, the disturbances will die out and finally be- 
come Mach waves. There may also be cases where a 
disturbance approaches a steady state, resulting in what 
is known as neutral stability. It is difficult to make 
these general statements more precise at this point; 
nevertheless, they can be accepted as qualitative defi- 
nitions of stability and instability. 

In the present investigation, the type of disturb- 
ances considered will be a finite linear pulse, moving up- 
stream, with a weak discontinuity at its head. By 
linear pulse is meant that at any instant the velocity- 
wave form is a linear function of space coordinates. 
This simplification may appear as artificial; however, 
observations in connection with ballistic waves show 
that the final state of any continuous pulse, traveling in 
still air, tends to become linear and is independent of 
the source of production. A continuous pulse prop- 
agating in a nonuniform steady flow may well develop 
into a complex shape, but a linear function could still 
be a good approximation to local conditions. 

If linear pulses are assumed, there will occur only 
two types: (1) a disturbance, propagating against a 
steady stream, whose ancestor was a continuous ex- 
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pansion pulse, will have a discontinuity at its head; 
(2) that disturbance whose ancestor was a continuous 
expansion pulse will have a discontinuity at its tail. 
These will be called compression and expansion pulses, 
respectively. Since an expansion pulse cannot pene- 
trate into the supercritical region but is dissipated at 
the sonic line, only the compression type of pulse will 
be considered. 

The problem of tracing the history of even a linear 
pulse is a formidable one. To reduce the amount of 
work involved, it is further assumed that all the dis- 
turbances were locally produced, by unsteadiness of the 
boundary layer, say, and are limited in extent later- 
ally. Thus, the triangular pulses encountered near the 
sonic line may be thought of as having originated in 
that neighborhood. Consequently, the motion of the 
pulse can be studied by solutions valid “in the small’’ 
instead of ‘‘in the large.’’ This brings about a great 
simplification because, by restricting the solutions to 
the neighborhood of the sonic line, one can introduce 
transonic approximations and obtain a simple solu- 
tion. 

Finally, it must also be said that, according to experi- 
mental observations in the case of airfoils, the pressure 
distribution over the surface does not exhibit the effects 
of separation or violent thickening of boundary layer 
until long after shock waves have been formed. This 
indicates that viscosity of the fluid has little influence on 
the character of the transonic flow outside the boundary 
layer, and therefore it can be neglected in the present 
study. Moreover, since the shock front is generally 
curved, the flow behind it cannot actually remain irro- 
tational and isentropic. However, by dealing with the 
initial stage of the motion, when the shock strength is 
weak and the vorticity generated is small; the flow can 
be considered as isentropic. If a disturbance increases 
in strength and tends to become a shock, the approxi- 
become worse and eventually break 


mation will 


down. 
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The future event inferred therefrom is tantamount to 
extrapolation. The validity of the conclusions drawn, 
of course, must be established by experiment. 


(2) DIFFERENTIAL EQUATION OF THE DISTURBANCES 


Let y’(X, Y, t) be the velocity potential of a two- 
dimensional motion, defined by 


, , 
u = gx’, VU = gy 
= a? — ao / (1) 
(u* + 0°) + - =o 
2 y-—I1 


Here, « and v are the velocity components, y is the 
ratio of specific heats, a is the local speed of sound, and 
dy is the value of a at stagnation conditions; the sub- 
scripts denote partial differentiation with respect to the 
coordinates X and Y and the time ¢. The differential 
equation satisfied by y’(X, Y, ¢) is 


§ € 9 9 , 
(a? — u*)p'xx + (a? — v*) o'yy — ¢'u — 


yy 
2uve’ xy = Que’ +, = 2Que’x, = O (2) 


If the motion is regarded as a superposition of a dis- 
turbance on a steady motion whose velocity potential is 
@(X, VY), the disturbance velocity potential @ can be 


defined as 

yg =G+¢ (3) 
For the present problem, ® is assumed known and satis- 
fies 
(c? — By*)Oyy + (C2? — Dy*)byy — 

2byPyPyy = 0 (4) 
with 
c? = ap? — [(y — 1)/2] (@x? + y’) 

By substitution of Eq. (3) in Eq. (2) and making use of 


Eq. (4), the differential equation governing the disturb- 
ance motion is found to be 


ct — ?,’ . (y a 1): _ (y = 1) (Pxodx + 2px") uot (y wo 1) (Po, 2 oy”) loxx + [c? ani ?,° ny 
(y - Ll): =e 1) (Pi dy = . op”) —tF + 1) (Pydy =i '/sy") |oyy = ee 2 (©, 9, + Pydy so 


?, oy os dy dy ldxy —* 2(Py, + by byt 


— 2b, + dx)oy: = 
(y — 1) (®yby + '/26)?) Onn + [Cy — Dede + ( 


[(y — ld + (y¥ + 1) (®xdx + '/ox") + 
— 1) (®xdx + . 2x") + ty + IP x 


(Py dy + ' apy”) |Pyy - 2[Px oy + Py dy + oxdy |®x, (9) 


For a given steady flow ®(X, VY), the problem would 
seem to be to solve this equation subject to a set of 
initial and boundary conditions. For an equation of 
such complexity, however, to seek an exact solution is 
not attractive; nor is it necessary because only certain 
features of the behavior of the disturbances are re- 
quired. For this reason, the problem will be solved by 
introducing appropriate approximations. 

It is assumed first that the disturbances are initially 
small and are produced in the neighborhood of the sonic 


line, so that the early stage of the pulse’s life is spent 
entirely in the transonic region. Asa result, the trans- 
onic characteristics of the steady flow will not be 
destroyed so long as the pulse is weak. Let € be the 
thickness parameter of the body. If the disturbances 
are small and their transonic character is preserved, 
then, according to transonic approximation,’ 0/OoX ~ 
0(1) and 0/OY & O(e”*). By the time the state of the 
pulse under consideration is reached, the wave front is 
practically straight, and its direction of motion is nearly 
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STABILITY OF 
parallel to the undisturbed flow. The time derivative 
must be the same order as the speed of propagation of 
plane waves. If the disturbances propagate against the 
steady flow, the time derivative must be of the order 
lby —c|. Now &x? — c? & O(e”’), therefore, 0/dt & 
Ole *), For the disturbances traveling downstream, 
the same argument shows that the time derivative 
must be of the order unity. 

As the pulse is predominantly plane, it might be rep- 
resented in the form g(x, y, rT) = g(x, 7) + eg) x 
(x, y, tT), Where the orders of magnitude of ¢ and ¢ 
must be alike in order to satisfy the boundary condi- 
tion of Eq. (9), below. This order of magnitude is 
arbitrary, depending on the strength of the disturb- 
ance considered. Since the motion must reduce to the 
propagation of a finite plane disturbance for vanishing 
¢, it will be assumed that ¢ and yg are O(e’’). In 
the following work, the two terms making up ¢(x, y, 7) 
will be handled together; nevertheless, the argument 
above serves to indicate the relative magnitudes of 
derivatives of g. To the order of e*, the upstream mov- 
ing disturbance satisfies 


(y + 1) (@e + Grrr + 262, + (CY H 1)Prr¢7 = 
tn ty 1 )$,,¢, — Grr (6) 
Similar consideration gives an equation 
Yrr = 26,2 = 0 (4) 
for the downstream moving pulses. Here, @ and ¢ are 
defined as follows: 


& = c*l(x + ) 
71) = c*ly 


and 


x= X/l, y= Y/l, r= ctl (8) 


where c* is the critical sound speed associated with the 
enthalpy of the steady flow and / is a characteristic 
dimension of the body. 

Eq. (7) shows that the downstream moving disturb- 
ances propagate without change of strength. If they 
are initially weak, they will stay weak. For this rea- 
son, only upstream moving disturbances require con- 
sideration. The equation governing the upstream mov- 
ing disturbances is a partial differential equation in three 
variables; aside from initial conditions to be specified 
later, the solution is required to satisfy the boundary 
condition: for y = Oandr >0O, 

Gy = $rE(x) (9) 
where g(x) is the slope of the body. It is assumed 
that the solution is regular on the axis and that on the 
body y is everywhere small. 


(3) WAVE ForM OF A PULSE 


Since, strictly speaking, Eq. (6) is valid only in the 
transonic zone, a local solution should be sufficient for 


TWO-DIMENSIONAL 


TRANSONIC FLOWS 3 


the present purpose. As the waves formed in the 
neighborhood of the sonic line must have lived long after 
the Riemann steepening has happened, the velocity 
wave form there can be approximated reasonably well 
by a discontinuous linear function in x and y at any 
time 7. If the origin of the axes coincides with the 
intersection of the sonic line with the body, the special 
wave form can be represented by 


l 
g = filr)x + fo(r)y + 5 [ fa(r)x? + 
2falr)xy + falr)y*] 


Moreover, 


(10) 


for values of x, y lying within the pulse area. 
®@,, too, is to be approximated by a linear function of 
x and y—namely, 


&, = ax + By (11) 
Of course, the sonic line, in general, is not a straight 
line; but, if the solution is restricted to a narrow region 
(a stream tube, say), Eq. (11) can certainly provide as 
accurate an approximation as desired by properly 
choosing the width of the stream tube. To be con- 
sistent with Eq. (11), the boundary condition (9) can 
be written as 

(12) 


oy = orm + kx), y = 0,7 >0 


where m and k are, respectively, the slope and curva- 
ture of the body at the sonic line. 

Upon substitution of Eq. (10), together with Eq. (11), 
in Eq. (6), there result 


(y + Ififs + 2f' + (vy + efi = fs 
(y + 1) (a + fs)fs + 2fs’ + (y + Deafs = 


—(y — I)afi’ — fr,’ (13) 
(y + 1) (8B + fadfs + 2fs’ + (vy + Doss = 
—(y — l)amf;’ — mf,” 
The same substitution in Eq. (12) gives 
he = mpi t (14) 


fa = kfi + mfs§ 


after the neglecting of terms of the order of x* and m?®. 
It is noticed that the order of the terms on the right- 
hand side of the second and third of Eqs. (13) is one 
order higher than that of the left-hand side. So, 
as a first approximation by neglecting high-order terms 
in f; and f,, the system possesses the following solu- 


tions: 
~ (Ae — 1)- Ki B ) —_ 
fi = (Ae” — 1) b ( —a)t+ k e | 
fo = mfi m 
fs = 2a(Ae”’ — 1)! (15) 
fs = (Ae”” — 1)! [28 + Be” *] 


fs = (2mv/k) [a — (8/m)\(Ae” — 1)! 


where v stands for (y + 1)a@ and A and B are constants 
of integration to be determined from initial condi- 
tion. It is to be noted that, in the case of one-dimen- 
sional motion, the solution does not reduce exactly to 
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that of Kantrowitz,'° inasmuch as fi(7) # 0 here in con- 
trast to the case considered by him. 

Having the velocity wave form chosen to be linear 
in x and y, it is natural to assume initially a pulse that 
is bounded at the front side by a discontinuity, in which 
the velocity decreases linearly from a given value to 
zero. Since the partial derivatives of x and y bear a 
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(12), the solution has to satisfy, at r = 0, 
¢r Sino — ¢,coso = a + bx (16) 


on the shock, where a is the shock angle and a, db are 
constants defining pressure ratio and its distribution, 
After approximating the shock front within a stream 
tube by a straight line perpendicular to the body, 





definite ratio, the rate of change in one direction fixes 
thatin the other. Asa result, the shape of the pulse can 
be determined by specifying the pressure discontinuity 
on the shock front or, alternatively, the normal ve- 
locity to the shock. Therefore, in addition to Eq. 


y— Yo = —(1/mMpo) (x — Xo) (17) 


where mp is the slope of the body at xo, a substitution of 
Eq. (10) in Eq. (16), making use of Eq. (17), yields two 
equations that possess a solution 


I 


9 3 / 
A 1+ — Jim — m*my — mM)a + = (a _— *) 


k(xo + myo) | m3 *) \ 
mok A \ k m is 1 + m? atid k (« om f 
2a 


- — ; moka + Bk (1 — mmo) + m? (« - *) vy— rp (18) 
m 


mok (1 + m?)'7A 
b ‘ 8 8 ] 
mom (1 + m*) la — — kB (xo + MYo) — m*(X + myo) (a — v |> 
a m ; m f 


9\ \/s 
(1 — mmo)a [ + m*-) 4 Xo + MYo ; 
9\'/2 rife 1? 
mo(1 + m?) k mo(1 + m?)”? 
Thus, for a given flow, the constants A and B may be regarded as functions of a and 6. By arbitrary choice of 


Being discontinuous in pressure, density, etc., the pulse cannot 
These provide an 


where 


A=- 


a and b, the motion within the pulse is defined. 
stay stationary but must propagate according to the Rankink Hugoniot shock-wave relations. 
It must be pointed out that, although the 
The discrepancies in pressure, 


additional equation to trace the successive locations of the —«1se. 
Rankine-Hugoniot relations are employed, the pulse is still regai 
density, etc., are permissible as long as the pulse is weak. 

Perhaps a word is also necessary concerning the process of continuation of solution to the next stream tube. 


d as isentropic. 


This is done in exactly the same manner except that the boundary condition at a stream line, say, is replaced by 
continuity in velocity. The continued solution, which possesses no special interesting feature, will be omitted 


here. 
(4) PROPAGATION OF A TRIANGULAR PULSE IN TRANSONIC REGION 


The solution can now be discussed without resort to detailed numerical calculation. Suppose the curvature k 


of the body is negative; the behavior of the pulse can be understood by examining the velocity components. 


& 2m (8 oie i 
or = (Ae” — 1)! ( —- a) + e’’’" + Yax + (28 + Be’ “)y 
1 
mB 


k \n k 
= (Ag” — 1)- Kk (° Z /) 4 “ ( 3 
oe k \m o k i k - m 4 


First, let (x, y) represent a decelerating flow so that the velocity gradients a and 8 are both negative. 
are to be discussed: a compression and an expansion pulse. For a compression pulse, the velocity slope 


y ») « 
g"'* 4+. (28 +. Be”’*) x 


6 


Two cases 
namely, 
Since x, y, 





[2a + m(26 + Be’’’*)]/(Ae’” — 1)—is positive and the velocity is negative—i.e., ¢, < O and gy, > 0. 
and m were assumed small, the character of the pulse depends upon the term {(2m/k) [(8/m) — a] + (B/k)e’"’*} X | 
(Ae”” — 1)~! and the slope 2a/(Ae’’ — 1). It follows that Ae’” — 1 must be negative and that (2m/k)[(8/m) — 
a] + (B/k)e’’” must be positive. 
Ae” — 1| decreases with time, while (2m/k) [(8/m) — a] + (B/k)e’”’? decreases or increases with time depending 
on whether B is positive or negative. In any case, since the rate of change of e’”’” is twice as slow as that of e”, 
the absolute value of the velocity is always increasing with time. 


Consequently, the momentum of the pulse 
continues to increase as the pulse propagates against the stream. 


Because of the fact that v is negative and 2m/k|(8/m) — a] is positive, | 
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It might be noted that the constant (2m/k) [(8/m) — a], being of the order €”’, is positive when m is negative. 


It might change sign if m were positive. 


This is not possible because the normal derivative of ®,, which is approxi- 


mately equal to 8 — am at the sonic line, is always negative in the case of a convex surface. 


Next, consider an expansion pulse. 
velocity, but there is a discontinuity at its ‘‘tail.”’ 
reverses its sign—that is, ¢, > 0 and gy, < 0. 
Asa result, B must be larger than 2(am — £). 
decrease with time and vanish when 7 is large. 


In this case, the ‘“‘head”’ of the pulse is continuous in both pressure and 
Here again, the slope of the pulse is positive, but the velocity 
This demands that (2m/k) [(8/m) — a] + (B/k)e’’’? be negative, 
This makes the absolute value of (2m/k) [(8/m) — a] + (B/k)e’"”? 
In other words, the slope of an expansion pulse whiile it propagates 


increases continuously to a steady value less than —2a, but its pressure discontinuity drops gradually to zero, 
Therefore, an expansion pulse, in the end, becomes a sound wave. 


Consider now the case of an accelerating transonic flow for which a is positive. 


then be written as 
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The velocity components can 
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Since v is positive, the expressions in the brackets will increase as e’’’“, and, as a consequence, the velocity will de- 


vr/2 


crease as ¢_ 
vious case. 
rate of its destruction. 


Furthermore, A — e~” tends to A instead of —1, this being an essential difference from the pre- 
Since A is a large number, it further reduces the final strength of the pulse and thus accelerates the 
In short, any pulse, regardless of whether it is compression or expansion, cannot be sus- 


tained long in an accelerating transonic flow once it is accidentally produced. 


(5) CONCLUSIONS AND DISCUSSION 

From the foregoing sections, certain conclusions can 
be drawn concerning the stability of mixed potential 
flows over convex surfaces. It is clear that the upper 
surface of a body is convex if k < 0. To be specific, 
the flow over the upper surface will be considered here. 
According to the results presented above, any dis- 
turbance, of the type considered here, produced in ac- 
celerating transonic flow over such a surface vanishes 
almost instantly, and the unsteady flow therefore ap- 
proaches the assumed steady flow as a final state. 
Thus, an accelerating flow of this type must be termed 
stable. 

On the other hand, a disturbance propagating against 
a decelerating flow grows continuously by consuming 
the kinetic energy of the steady flow. 
of converting kinetic energy into pressure will continue 


This process 


to prevail until the pulse is so strong that it changes 
drastically the character of the whole flow pattern. 
The development of such a pulse is illustrated in Fig. 1. 
The decelerating transonic flow, therefore, is unstable 
to compression pulses. From this, one concludes that 
smooth transonic flows with local supersonic regions 
involving both accelerating and decelerating flows are 
unstable. 

It is interesting to see that the curvature of the body 
plays no decisive role in the question of stability of two- 
dimensional smooth mixed flows. The major effect is 
evidently that of the velocity gradient ®,,, as in one- 
dimensional flow. One may suspect that the constant 
(2m/k)((8/m) — a] might have an important influ- 


ence. In the neighborhood of the sonic line, it can be 
shown that it represents the ratio of the normal deriv- 
ative of ®, to the curvature %. As long as the slope of 
the surface there is small, the normal derivative of %, 
By the 
irrotationality condition, it follows that this constant 


is approximately the same as that of the speed. 


There- 
fore, if the curvature is varied keeping other conditions 


is practically independent of the curvature k. 


unchanged, the constant will remain the same, and, as 
a result, the pulse will not differ. 

As stated above, the major effect that decides the 
question of stability or instability is the velocity gra- 
dient. The conclusion drawn, therefore, is valid only 
Presumably, this 
But the 
Is it possible that &,, does vanish and 
that , has a point of inflection there 
an*+...? 


is still open, and the problem has to be re-examined. 


if ®,, is finite and nonvanishing. 
covers the most important cases in practice. 
question arises: 
namely, ®, = 
If it does happen, the question of stability 


A word is also necessary regarding the thickness effect 
of the body. For convenience, the approximation is 
based on a small thickness parameter. 
proximation is applied only locally, the solution will be 
valid as long as the slope of the surface in the region 


Since the ap- 


considered is small, regardless of whether the body is+ 
thick or thin. In the case of a thick body, the only 
restriction would be for low free-stream Mach Num- 
ber, so that the local supersonic region occupies the flat 
portion of the surface. Therefore, basically, in con- 
sideration of the problem of stability, the thickness 
effect, too, does not enter. 
(Continued on page 54) 
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The Characteristics of Supersonic Wings 
Having Biconvex Sections 


BEVERLY BEANE* 
United Aircraft Corporation 


SUMMARY 


The wave drags of a representative series of thin, tapered 
plan-form wings having multislope surfaces approximating the 
circular are or biconvex section are computed by means of a 
generalization of the linear source-sink solution of reference 1. 
The computational methods are described, and the results are 
compared with those of similar wings having double-wedge sec- 
tion. Limited experimental checks of the computed resuits are 
included. It is concluded that the linear theory gives a reason- 
able first-order approximation of the characteristics of biconvex 
section wings and that for most applications the biconvex section 
will be superior to the double-wedge section. 


INTRODUCTION 


- 1946, PUCKETT PRESENTED a solution for the wave 
drag of thin delta plan-form wings of double-wedge 
section.'. This paper was followed in rapid succession 
by solutions for the lift of such wings® * and for the lift 
and drag of other plan forms of practical interest.‘~*! 
These solutions, most of which are well summarized in 
22, completely describe the major char- 


wey 


reference 
acteristics of supersonic wings having double-wedge sec- 
tions. 

This paper is concerned with the characteristics of 
supersonic wings having biconvex or circular-are sec- 
tions. 

All of the finite-wing solutions are of the so-called 
linear form and are based on the assumptions that the 
fluid is a homogeneous, nonviscous medium conform- 
ing to the ideal gas laws, that the disturbance pertur- 
bations are small enough so that their squares and 
products can be neglected relative to the square of the 
stream velocity, and that the boundary conditions im- 
posed by the wing may be satisfied in the horizontal 
plane rather than on the actual airfoil surface. On the 
basis of these assumptions, the potential equation may 
be linearized and solved by any one of several available 
methods.!: Of these, the source distribution 
method is the most widely used and is the one herein 
employed. In this method the thickness profile of 
the wing is simulated by a source distribution in the 
the angle of attack, by a similar 
The pressure and forces at each 


23 —29 


horizontal plane; 
doublet distribution. 
point on the wing are obtained by integration of the 
effect of these source and doublet distributions. 


Presented at the Aerodynamics Session, Eighteenth Annual 
Meeting, I.A.S., New York, January 23-26, 1950. 
* Analytical Engineer, Research Department. 


Most wings satisfy the assumptions of the linear 
theory and experiment has, to a large degree, verified 
the validity of its overall results. 

Fig. 1 presents a reference summary of the finite- 
wing solutions. As may be noted, while the character- 
istics of wings of double-wedge section have been ex- 
tensively studied, comparatively little has been done 
on the characteristics of wings of biconvex section. 
With the exceptions of Jones’ solution for the drag of 
sweptback constant chord wings of biconvex section'! 
and Puckett’s investigation of the effect of position of 
maximum thickness, ! all of the wave drag solutions have 
been for wings of symmetrical double-wedge section. 
Unfortunately, this section has poor structural char- 
acteristics, so that these solutions do not correspond to 
the sections in which the designer is most interested. 

The double-wedge section was chosen for these initial 
solutions primarily because of the resulting simplifi- 
cation of the solution and not because of any inherent 
advantage of the profile. Two-dimensional section 
characteristics alone strongly suggest that, because of 
its better structural shape, the biconvex or circular- 
are section would be better for a practical wing de- 
sign (see Fig. 2). With the exception of the compara- 
tively simple case of the constant chord wing, no solu- 
tion has yet been published for finite wings having bi- 
convex or any smoothly curved profile. This is be- 
cause of the difficulty that arises in the solution of the 
resulting integrals. Thus, although on the basis of two- 
dimensional section data it has generally been felt that 
the biconvex section would probably be advantageous, 


















































SUMMARY_OF FINITE WING SOLUTIONS 
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SUMMARY __OF TWO _ DIMENSIONAL 
SECTION CHARACTERISTICS 
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a direct comparison could not be made because drag 
solutions fgr wings of biconvex section were not avail- 
able. 

A curious situation has thus arisen. Theoretical 
analyses have emphasized the double-wedge section 
because of its simplicity. 
menters, in turn, have concentrated on it because of 
The result is 


computational Experi- 
the existence of the theoretical solutions. 
that today we have extensive theoretical and experi- 
mental drag data on a section that is probably of second- 
ary interest. 

It has always been recognized that the drag of 
smoothly curved profiles could be computed to any 
desired accuracy by approximating the section by a 
many-sided polygon, which, in turn, may be analyzed 
using a generalization of the solutions developed in 
reference 1. Through an inelegant and_ laborious 
method, the desired solution can be thus obtained. 
While it is to be hoped that more concise treatments 
will eventually be worked out, for the present this 
method serves at least an immediate purpose. A 
sufficient number of such solutions have 
out so that it is possible to present a fairly complete 
picture of the characteristics of biconvex wings. While 
it is not maintained that this section is optimum, this 


been carried 
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work substantiates that, for most applications, the bj- 
convex section is definitely superior to the symmetrical 
double wedge. 

It is the purpose of this paper to present these re. 
sults and to compare them with the previously pub- 
lished results for wings of double-wedge section. 

The type of plan form and range of variables for which 
solutions have been worked out are shown in Fig. 3. 
These include the symmetrical diamond, delta, and ar- 
row plan forms. Since most of the variables are noted 
in nondimensional form, the range of leading-edge sweep 
angles at a representative Mach Number of 2 are indi- 
The First Reverse Flow Theo- 
rem”*~* states that the results of these analyses are 


cated in parentheses. 


equally applicable to geometrically similar swept- 
forward wings. 

The method used in obtaining the drag of these wings 
is described in the next section. Although here applied 
to investigate the characteristics of biconvex section 
wings, the method as outlined is applicable to the pre- 
diction of the drag characteristics of thin pointed-nose 


sections of arbitrary shape. 


NOTATION 


A = area 
k 
A s..> 
|? 
ky, ' : 
PD ny, k = functions representing pressure integrals 
j ky, 
Fi Nas 
7 


B = cotangent of Mach angle 

Cp = drag coefficient 

C, = pressure coefficient 

oa = lift-curve slope 

da 

Cr = root chord 

Cp = Projected chord [ce/cg = 1 — (kw +1/hki)] 
y = leading-edge sweepback angle 
C; = tip chord 

k = tangent of sweepback angle 

ky = tan, 

Rn +1 = tang 





RANGE OF VARIABLES FOR WHICH SOLUTION CARRIED OUT 
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SUPERSONIC WINGS 

r. = surface slope (relative to direction of flow) 

d = taper ratio (A = ¢;/cr) 

My = free-stream Mach Number 

m = tan '8 (complement of Mach angle) 

n = sweepback parameter (n = k/B) 

N = number of slopes into which plan form is 
divided 

qo = dynamic pressure in free stream [go = 
(1/2) po Uo? | 

qs = source strength (g, = Uod,/7m) 

5 = kn + i/ki 

o = trailing-edge sweepback angle 

tg = integration variables 

u = horizontal disturbance velocity 

Up = free-stream velocity 

Subscripts: 

B = biconvex 

W = double wedge 

© = infinite aspect ratio 


METHODS OF ANALYSIS 


The fundamental linearized equations for deriving 
the pressure drag of a thin airfoil at supersonic speeds 
are presented in reference 1. The method consists 
of solving for the disturbance velocity potential due to a 
sheet of sources distributed over the area of the wing, 
imposing the boundary condition that the flow must be 
parallel to the airfoil at its surface. The pressure at 
any point on the wing is then proportional to the flow- 
direction derivative of the velocity potential. 

A similar procedure would be the most direct ap- 
proach to the problem of finding the drag of thin bicon- 
vex wings, but the continuous curvature of the airfoil 
surface introduces an additional factor that makes the 
integration of the potential equations extremely diffi- 
cult. For this reason, the present paper considers a 
different, but simple, approach, approximating the 
biconvex profile (or any curved profile) by a many- 
sided polygon and obtaining drag-coefficient solutions 
by superimposing the known solutions for single trian- 
gular wings of constant slope. This is the identical 
procedure followed by Puckett and Stewart in order 
to calculate the drag of two- and three-slope wings. 

Pressure distributions over single delta wings of 
constant surface slope A, are derived in reference 1, the 
boundary conditions being satisfied by a distribution of 
constant strength sources within the delta. These 
solutions are applicable to any triangular wing, pro- 
vided the trailing edge Mach cone does not intersect 
the wing. Since solutions to linear equations may be 
added, any arbitrary distribution of slopes over an air- 
foil may be obtained by superimposing several such 
triangles with common trailing edge (see Fig. 4). The 
resultant pressure distribution is then made up of super- 
posed pressure distributions, each associated with one 
of the triangles having vertices at A, B, C,..., ete., 
and having a common trailing edge at OA’. 

It is readily seen that a biconvex airfoil can be ap- 
proximated by a proper superposition of triangles of 
appropriate source strength and could be simulated 
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Fic. 4. Multisloped triangular wing plan form. 
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Single-slope triangle. 
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exactly by superimposing an infinite number of such 
triangles—i.e., by allowing the lengths AB, BC, CD,..., 
etc. (Fig. +), to approach zero. 

Considering a representative single-slope wing such 
as Fig. 5, the linear theory indicates the pressure is con- 
stant along rays through the vertex A; that is, along 
lines ¢ = constant, where t = ky/x and k is the tangent 
of the sweepback angle. Defining the pressure co- 
efficient by C, = —(2/Uo)u, where Up = free stream 
velocity and u = horizontal disturbance velocity, 
the following relations are derived in reference 1: 

(A) For pressures within the triangle when the Mach 
wave is ahead of the leading edge (i.e., n > 1 and 
F< 2), 

tr, n?* — -? 


cosh! (1) 


C, = 
: | 1 — /? 


BrvV n? — 
(B) For pressures ahead of the triangle but behind 
the Mach wave when the Mach wave is ahead of the 


leading edge (i.e.,n >¢> 1), 


n> — | 
cosh~! \". (2) 
2 


(C) .For pressures within the triangle and behind the 
Mach wave when the leading edge is ahead of the Mach 
wave (i.e.,4< n< 1), 


4X, 


i = Fr 
Brv n° — 
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The above basic pressure coefficients multiplied by 
the local slope of the airfoil and integrated over both 
airfoil surfaces will now yield the drag, or 

dD/qo = C,d,;'dA (5 
where go is the dynamic pressure in the free stream 
[go = (1/2)poUo"| and X,’ here is the local slope of the 


airfoil at the element of area dA. Consider now a 
multiple-slope wing as in Fig. 4, and let A; be the slope 





in the region A, (triangle OAB); 2, the slope in the 


| —— X region A» (triangle OBC), ete. It is convenient, for 
AC, integration purposes, to consider the wing as _ being 


composed of several superposed source distributions, 





Fic. 6. Elementary integration area (f < 1). 
each constant within a complete triangle back to the 


trailing edge OA’. There is then one pressure distribu- 





tion due to the sources of strength gq, = UodA:1/7, dis- 
tributed over area OA’A; a second pressure distribu- 
tion due to the sources of strength g. = Uo(A2 — Ai) /7, 
distributed over area OA’B; a third pressure distribu- 
tion due to sources of strength gz; = Uo(As — As)/7, 
distributed over area OA’C; etc. In computing the 





net drag, it is necessary to take into account all of these 
superimposed pressure distributions and therefore to 
remember, in particular, that A,’ in the expression for 





dD/qo is the actual local slope of the elementary area 


x dA whose drag contribution is being computed, while 
|— 16, —-| A, in the C,-equations is the slope associated with the 


source strength of the particular triangle whose effect 





Fic. 7. Elementary integration area (f > 1). ‘ : : é i ; 
is being added. (For a single-slope wing, A, and \, 


are, of course, equal.) 


tr ee Two elementary area expressions are needed: 

,= , 5 ~ sin ; . (3) (A) For the integration of the basic pressure dis- 
/ ad ) P- — t? ‘ ° om ‘ 5 

BrvV | n tributions of Eqs. (1), (3), and (4), over areas of the type 


(D) For pressures within the triangle and ahead shown in Fig. 6, the elementary area is 


of the Mach wave when the leading edge is ahead of 1 Ace? dt, 
the Mach wave (i.e.,7 <t< 1), dA = 2 ky [l — (By/Ra)tal? (6) 
C, = 2d,/BV 1 — n? (4) (B) For integrating the pressure distribution of 


; / Mt ; Eq. (2), over an area such as Fig. 7, 
In these equations B = M,? — 1 (Mp is the free- 


stream Mach Number), 1 = k/®8, and X, is the surface talk ses 1 Acr? dt, (7) 
slope of the wing. ~~ 2 hy [(Ra/Re)ty — 12 ' 
The required integrals are: 
(A) For the basic pressure of Eq. (1) in an area of the type shown in Fig. 6, 
2r, Acr* l , Na” — ty? dtq 
| C,dA = = { cosh~! ,/| = - 
Ag Br Rag Vn. —1JS0 V1l-4@? [1 — (Rp/Ra)tal ? 
= (2A,/Bm) (Acy?/Ra) B(na, Ry/Ra) (8) 
where 
ky l In 1, k, cosh! n, 2 V1 — ny? 
B (ma ) = a | + + tan—-! — 3" .  < I 
a 1 — (ky/ka)? Ln? — 1 ka Vn — 1 V/1 — ny? Na — NM + VA = | ' 
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or 


P ( *) 1 | In mq 4 k, cosh! ny 4 
Nay - 9 > ° 
ka 1 — (Rp ke)? Vn, — | Ra Vn? — ] 


I 2 ny? — |] 
; In (: + )} nN, > 1) Ya) 
Vn? — 1 Ne — ty + Vn? — 1 — Vn? — 1 Me 


(B) For the basic pressure distribution of Eq. (2) in an area of the type of Fig. 7, 


= Zr ACr’ l "b Ny? — 1 dt, 
C, dA = r cosh~! ; Z 
J A, Br ky Vv Ny? —l1/J/1 ty? — | [(Ra/Rr)tr = 1} ? 


= (2,/Ba) (Ace?/Ry) (Ro/Ra)?D(mo, Ry/Ra) (10) 
where 


o( f l i Inm l ; |" — (Ry/Ra) + V [my? — (ky/Raq)?| (mM? — | (11) 
Nn, = ; n 
we 1 — (ky/Ra)? \\/my?—1 Vn? —1 m[1 — (Ro/Ra)] f 


(C) For the basic pressure distributions of Eqs. (3) and (4) in an area of the type of Fig. 6, 


2r, Acr? 1 ; ity m . a? — tq? Ut, 
I CdA = — - ; [ id —— f sin! ‘ — A 
A, Br ka V1 —n2 \Jo 2 [1 — (ko/Radtel? S02 V1 — 14,2 [1 — (be/Ra)tal ?9 


= (2d;/Bx) (Acp?/Ra)f(me, Ry/Ra) (12) 


where 


i ( *) : | — ‘nt ( + sin )| (13) 
Na; = cos~! mq sin~! n : 
ka) 1 — (h/ka)? LW — mg? V1 — m? \2 ; 


It is now possible to compute the drag of the multislope wing of Fig. 4 by using appropriate combinations of 


these functions. The total drag will consist of contributions in regions A), A», A3,..., and Ay due to the source 
distributions in triangles OA'A, OA'B, OA'C,..., ete. 
Let 


Dy = drag in region A, 
Dy = drag in region A» 


A 
due to sources gq; = nin OA'A 
T 
Diy = drag in region A y J 
Dy = drag in region A, ) 
| 
Uo(A2 — Ai). 
» due to sources g. = s - in OA'B 
| T 
Dox = drag in region A y } 
Dy, = drag in region A; ) 
Dy2 = drag in region A» | 
| Uo(Aw — Aw-i) . 
‘ due to sources gy = alal v—" in OA'B’ 
TT 


region Ay 


_e 
~ 


— drag 11 
Then 
Drota = (Du + Dis + Dis +... + Din) + (Dn + Dn + Dg +... + Day) + ~~ + 

(Dm + Dye + Dnz3+...+ Dnwn) (14) 


The form that these drag components take depends, of course, upon the position of the Mach waves relative 
to the leading and trailing edges of the region whose contribution is being computed—that is, upon the value of 
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the ratios k,/8, ke, etc. For the general case, ki > ky >k3 >... > Rp > B > Roy >... > Ry4i—that is, all the 
ridge lines up to (and including) the pth are swept at an angle greater than the Mach angle, while the remaining 
ridge lines are swept less than the Mach angle (/ is equal to zero when the wing leading edge is ahead of the Mach 
cone). The drag components are, then, 
Dy 4i,? AB? ( =) 
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In summing these dtag contributions, it is convenient to consider that the regions of various slope are all equal 
= Ay,andsoAB = BC = CD =... = XpXp41 = --- = BA’ = 


2)crCp/k,| is, then, 


in area—that is, A; = Az =...=A,=... 
ce N. The total drag coefficient based on actual wing plan area [Ar = (1 
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In using this equation to compute the drag of a multislope wing, it must be remembered that the derivation in- 
cludes the assumption that the ridge lines all intersect the base chord at equal intervals (Acg = ce’ N). For cases 
where this assumption is not true, it is a relatively simple matter to obtain applicable equations by adding the 
components in Eq. (15), using the appropriate values of AB, BC, CD, etc. However, it is possible, and probably 
(16) directly, considering the plan form divided by imaginary ridge lines into regions of equal 


simpler, to use Eq. ( 
Terms involving the & ratios of the 


area with the correct slopes and substituting the corresponding value of V. 
imaginary ridge lines will then go out, since the \ differences are zero. 
Setting V equal to infinity in Eq. (16) provides an exact solution to the problem of the drag coefficient of a bi- 


convex airfoil, but attempts to sum the series represented were unsuccessful. However, finite values of V may 








14 JOURNAL OF THE AERONAUTICAL SCIENCES—JANUARY, 1951 


be used to calculate approximate solutions whose accuracy can be increased by increasing N as far as computa- 
tional time permits. 

For small thickness ratios, the biconvex airfoil can be accurately represented by the symmetrical biparabolic 
airfoil. If the section composed of two parabolic arcs is then assumed approximated by an inscribed polygon 
with 2 faces of equal horizontal projection, the geometric relationships 


A = —Ay = 2t/o)[(N — 1)/N], A— Ay = AA = (4/N)(t/c) 


can be used to reduce Eq. (16) to the following forms: 
(A) When the initial Mach wave is ahead of the wing leading edge, 
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(B) When the leading edge lies ahead of its initial Mach wave, 
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It is the conclusion of this study that it is probably considered to be divided into 16 slopes in order to 
not worth while to use an N value greater than about 16 keep the number of calculations within reason and yet 
(see Fig. 8 for illustration). to approach satisfactorily the true drag value for a 

Thus, to compute the wave drag of the representa- smooth are profile. For the symmetrical biconvex 
tive plan forms included in this study, the profile was profile, N [in Eq. (17)] was then equal to 16. To 
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proximated by symmetrical 2.\-sided polygon). 


investigate the effects of maximum thickness position, 
profiles were selected which were composed by two 
circular semiarcs different radii, tangent (hori- 
ontally) at the point of maximum thickness. To 
facilitate calculation, the profile was again divided 
into 16 slopes approximating the true arcs, eight 
forward of the maximum thickness position and eight 


of 


alt. 

As stated previously, the pressure drag of a wing 
is proportional to the flow-direction derivative of the 
velocity potential due to a sheet of sources distributed 
over the area of the wing. The pressure at any point 
P (at zero angle of attack), for instance, is determined 
only by the sources within an area bounded by the 
Mach lines forward from P and the wing boundaries. 
The effect of a finite tip is thus confined to the region of 
the wing outboard of the Mach wave springing from 
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= — cosh 
BrvVn,” — | 


[Eq. (18S) assumes that the Mach waves do not cross 
on the surface of the wing.| The wave drag contribu- 
tion of the tip is then the integral of Eq. (18) multiplied 
by the local slope of the airfoil and integrated over the 
area enclosed by the tip Mach For a multislope 
airfoil, the pressure distributions can be superimposed 
It is found that the net tip drag of a sym- 
airfoil is zero. 


yave. 


as before. 
metrical biconvex or double-wedge 
Thus, an equation for the drag of tapered plan forms 
with biconvex sections can be set up similar in form to 
Eq. (17), but the pressure functions now represent the 
integrals of the pressure coefficients of Eqs. (1) through 
(4) over tapered areas, such as Fig. 9. 


RESULTS 


The final values of Cp[6/(t/c)*| for a symmetrical 
biconvex airfoil, as given by Eq. (17), are functions of 
three parameters: n, the leading-edge sweep ratio; s, the 


trailing-edge sweep ratio relative to that of the lead- 





HAVING 


1 


BICONVEX SECTIONS 





/ 


x, =— / 
Area of 
missing sources -+-~ nN / me 





x' 











Le Am | 
r ACR mn 
Fic. 9. Area affected by tip cutoff 


the forward edge of the wing tip. The pressure distri- 
bution over any swept tapered wing will therefore be the 
same as that on the corresponding triangular wing 
(taper ratio zero), except that, over the part of the wing 
enclosed by the tip Mach wave, the pressure increment 
due to the ‘“‘missing’’ source elements (see Fig. 9) 
must be evaluated and subtracted from the distribu- 
tion given by Eqs. (1) through (4). 

The pressure increment due to the missing sources 


is found to be constant along rays through the wing 


tip—.e., along lines g = constant where g = By’/x’ 
(Fig. 9). The resulting tip-effect pressure coefficient 
is 
; Na + g 
sin! & > Mock B) 
+ Mag 
(18) 
Na + g 
» (t,. > 1) 
1 + Mag 


ing edge; and X, the taper ratio. In order to present 
a fairly complete picture of the characteristics of bi- 
convex section wings, wave drags have been computed 
for three representative wing plan forms (Figs. 10-12). 
Lift-curve slopes, although independent of cross section, 


are included in the figures for completeness. 


Fig. 10 presents the results for symmetrically tapered 
The drags are presented in terms of C,[8 + 
which physically represents the ratio of the 
wave drag coefficient to that of a two-dimensional 
double-wedge section of equal thickness ratio. Re- 
sults are plotted against n, the ratio of the tangent of 
the leading edge sweep angle (vy) to that of the leading 
This type of presentation is 
independent of thickness ratio Mach Number. 
Results are shown for taper ratios (i.e., tip chord to 
root chord) 0.0, 0.25, and 0.50. Corresponding 
double-wedge drags* * * are included for comparison. 
It is seen that the drag increases as the sweep angle 
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Fic. 10. 
is increased and that the effect of cutting off the tip 
is slight. In all cases, for equal thickness ratio, the 
drag of the biconvex section wing is greater than that 
of the double-wedge section wing. The comparison 
between the sections will be considered in greater de- 
tail later. 

The lift-curve slopes for the triangular plan forms 
were obtained from reference 4. For taper ratios other 
than zero, the slopes were computed using (1) the re- 
sults of reference 15 for the case of wings with sub- 
sonic leading edge (a close approximation found by 
application of the method of references 31 and 32), and 
(2) an original analysis based upon the method of 
Evvard*' for the case of wings with supersonic leading 
edge. In all cases, both the lift and the wave drag 
solutions were limited to wings that have supersonic 
trailing edges and for which the Mach lines emanating 
from the wing tips do not intersect on the wing. 

Tables 1 and 2 were prepared so that the results 
plotted in terms of C,[6/4(t/c)*| and (dC,/da) (8/4) 
may be readily converted into actual coefficients of wave 
drag and lift for a given plan form and Mach Number. 

Fig. 11 presents the results for wings of delta plan 
form. Again the drag is seen to increase with sweep 
angle until the Mach line crosses the leading edge 
1.0), after which the drag falls with further 
The biconvex wing again follows, 


G2..sts = 
increases in sweep. 
on a higher level, the same general trend as the double- 
wedge wing. 

Finally, in Fig. 12 are shown the results for a typical 
arrow plan form. The trends are everywhere the same 
as for the other plan forms. 

In all cases it must be remembered that the drag 
coefficients shown are wave drag and not total drag 
and that the lift-curve slopes are isolated wing char- 
acteristics and do not include body upwash effects. 

Fig. 13 shows the effect of varying the position of 
maximum thickness of a representative triangular 
plan form. It may be noted that, for minimum wave 
drag of wings swept well behind the Mach line, the 
maximum thickness of the circular-are profile should 
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theoretically be located well forward of the mid-chord 
position. This same conclusion was also indicated by 
theory for plan forms with double-wedge cross section. ! 
For that case, however, experiment® has shown that 
the sharp ridge line, when located well forward, causes 
an earlier transition of the flow and increases the skin 
friction enough to offset the predicted decrease in wave 
drag. It might be supposed, however, that, in the case 
of the wing with circular-are profile where there is a 
gradual change in curvature rather than the sharp ridge 
line of the double wedge, this forward flow transition 
will not occur and a decrease in total drag may actually 
be obtained. 

The wave drag characteristics of plan forms with 
biconvex and with double-wedge sections are directly 
compared in Fig. 14. As expected, this figure shows 
that, on the basis of equal thickness ratios, the drag of 
a plan form with biconvex section is higher than that 
of the same plan form with double-wedge section for 
all sweepback angles. However, comparing the two 
sections on the basis of equal moment of inertia about 
the structural criterion of equal 
the biconvex wing has a lower 


the chordwise axis 
bending deflection 
drag. Comparisons on the basis of equal section 
modulus, equal area, or equal moment of inertia about 
) axis lead to still greater advantages for 
It is on the basis of these results 


the major han 
the biconvex wing. 
that it is concluded that the biconvex section will, for 
most applications, be superior to the double-wedge 


section. 


COMPARISON WITH EXPERIMENT 


In order to check the validity of the conclusions 
drawn from the theoretical analysis presented herein, 
wings were tested in the 
The test setup is shown in 


a limited number of 
U.A.C. supersonic tunnel. 
Fig. 15. Sample results from these tests are shown in 
Fig. 16. 

Fig. 16(a) shows the drag of an arrow-type biconvex 
wing with maximum thickness at 50 per cent chord. 
The solid curve is the theoretical calculated wave drag, 
including angle-of-attack drag (i.e., Cha). The dashed 
curve is the calculated wave drag plus laminar skin 
friction (Cp, = 0.0053), and the dash-dot curve is wave 
drag plus turbulent skin friction (Cp, = 0.0096). It 1s 
seen in this case that the computed drag at zero angle 
of attack is in satisfactory agreement with the test 
point, and the indication is that the boundary layer 1s 
largely laminar. As shown later, all of the tests did not 
show such good agreement. 

Fig. 16(b) shows the result for a wing of geometrically 
identical plan form and thickness ratio to that shown 
in Fig. 16(a) but having its maximum thickness at 
20 per cent chord. This wing was tested to determine 
whether the theoretically predicted drag reduction 


could be realized. It is seen that the predicted drag 
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reduction is realized, the proportion of skin friction 
drag appearing to be about the same as for the wing 
It thus 
appears, on the basis of this limited check, that the use 


with maximum thickness at 50 per cent chord. 


of the biconvex section, eliminating the double-wedge 
ridge line and its associated pressure discontinuity, 
may make possible the attainment of the theoretically 
predicted drag reduction associated with forward loca 
tions of the maximum thickness line. 

Fig. 16(c) shows the result for another arrow-type 
biconvex section wing with maximum thickness at 
90 per cent chord. For this test the Mach line lay 
along the leading edge of the wing. It may be recalled 
that, for this condition, the theory predicted a drag 


peak. Itis shown that this drag peak is not realized in 
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typical arrow plan form. 
practice. The drag is, in fact, of the same order of 
magnitude as for cases in Figs. 16(a) and 16(b). 

In general, all of the tests showed drags in agreement 
with, or less than predicted by, the theory. Fig. 
16(c) was the poorest agreement obtained. 

All of the experimental results are summarized on 
Fig. 14, where they are superimposed on the theoretical 
results for comparison. Allowing for the effect of 
the linearizing assumptions on the theoretical results 
it appears from these limited checks that the linear 
theory gives a reasonable approximation of the char- 
acteristics of biconvex section wings. It also appears 
that the conclusion as to the superiority of the bicon 
vex section over the double wedge will be substantiated 


in practice. 
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The Rolling Up of the Trailing Vortex Sheet 
and Its Effect on the Downwash 


Behind Wings 


JOHN R. SPREITER* ann ALVIN H. SACKS* 
Ames Aeronautical Laboratory, N.A.C.A. 


SUMMARY 


The motion of the trailing vortices associated with a lifting wing 
is investigated by theoretical and visual-flow methods for the 
purpose of determining the proper vortex distribution to be used 
for downwash calculations Both subsonic and supersonic speeds 
are considered in the analysis. 

It is found that the degree to which the vortices are rolled up 
depends upon the distance behind the wing and upon the lift 
coefficient, span loading, and aspect ratio of the wing. While the 
rolling up of the trailing vortices associated with high-aspect 
ratio wings is of little practical importance, it is shown that, with 
low-aspect-ratio wings, the trailing vortex sheet may become es 
sentially rolled up into two trailing vortex cores within a chord 
length of the trailing edge 

The downwash fields associated with the two limiting cases of 
the flat vortex sheet and the fully rolled-up vortices are investi 
gated in detail for both subsonic and supersonic speeds. The 
intermediate case in which the rolling-up process is only partially 


completed at the tail position is also discussed 


INTRODUCTION 
ee egaeonen STUDIES OF WINGS of finite span at 
both subsonic and supersonic speeds are generally 
based on the assumption that the perturbation veloci- 
ties are much smaller than the free-stream velocity. 
The differential equation for the perturbation velocity 
potential y is then 


(1 — M?) gee + Gy + Oe = 0 (1) 


where .\/ represents the free-stream Mach Number and 


the coordinate system is as shown in Fig. 1. In cer- 
tain important cases, however, (1 — J/*)gy, is much 


smaller than ¢,, and ¢., and Eq. (1) may be reduced to 
the form 


Puy = ¢2 = 0 (2) 


With low-aspect-ratio wings, conditions permitting the 
use of Eq. (2) are found in the vicinity of both the wing 
and the wake at all Mach With high- 
aspect-ratio wings, Eq. (2) is valid at large distances 
behind the wing but may not be used generally in the 
Another important case for 


Numbers. ! 


vicinity of the wing. 
which Eq. (2) may be used is that for the flight of swept 
lifting surfaces at sonic velocity, since, then, the first 


Condensed version of a paper presented at the Aerodynamics 
Session, Eighteenth Annual Meeting, I.A.S., New York, January 
23-26, 1950. Received July 28, 1950. 

*Acronautical Research Scientist 


term of Eq. (1) is omitted because the coefficient 


(1 — AJ?) 1s zero.” 4 


The boundary conditions are prescribed on the wing, 
usually assumed to lie in the xy plane, either by speci- 
fying ¢. by giving the shape of the wing surface or by 
specifying yg, by giving the load distribution. In addi- 
wake. At this point, it has become customary at both 
subsonic and supersonic speeds to admit the assump- 
tions originally introduced by Prandtl* ° in the estab- 
lishment of lifting-line theory. They are (1) that the 
vortex wake of finite thickness may be replaced by an 
infinitesimally thin surface of discontinuity, designated 
the trailing vortex sheet, and (2) that the trailing vor- 
tex sheet remains flat and extends downstream from the 
wing in the free-stream direction. With these assump- 
tions, the boundary condition in the wake may be ex- 
pressed as a discontinuity in the potential Ags having 
a spanwise distribution identical with that at the trail- 
ing edge of the wing. 


Although it has been firmly established that these 
assumptions are sufficiently valid for the prediction of 
the forces and moments on finite-span wings, the ex- 
tent of their validity for the calculation of the downwash 
Prandtl 


pointed out that the vortex sheet begins as a flat sheet 


behind wings has never been clearly defined. 


at the trailing edge and rolls up into two vortex cores at 
great distances behind the wing. Thus he, and later 
Glauert,® advocated the use of a flat vortex sheet for 
calculating the downwash for points near the wing and a 
single horseshoe vortex for points far behind the wing. 
Helmbold’ * recognized the influence of the lift co- 
efficient and recommended the use of the flat vortex 
sheet for lift and the horseshoe 
vortex for large lift coefficients. Since later experi- 
mental studies,* ' together with Kaden’s theoretical 


small coefficients 


analysis,'' showed that the vortex sheet was nearly flat 
at the tail position throughout almost all of the useful 
lift-coefficient range of the unswept high-aspect-ratio 
wings in use at the time, most subsequent theories 
adopted the assumption that the vortices were arranged 
in a flat sheet. Since the early works were almost in- 
variably confined to unswept wings of high aspect 
ratio, an important gap in our knowledge exists for the 
low-aspect-ratio wings and swept wings being intro 
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View of wing and coordinate axes. 


duced for high-speed flight. This study will attempt 


to supply an understanding of some of these problems. 


A fundamental difficulty in the full analysis of down- 
wash lies in the fact that the paths of the trailing vor- 
tices cannot be prescribed a priort, as is done in classical 
wing theory, but must coincide with the stream lines, 
which are in turn influenced by the vortices. The 
present treatment will serve as an approximation in 
that a knowledge of the shape of the vortex sheet will 
be obtained through some two-dimensional, as well as 
three-dimensional, considerations; the downwash will 
then be determined by means of the customary three- 
dimensional analysis. 


SHAPE OF VORTEX SHEET 


Rate of Rolling Up of Vortex Sheet 


Considerable insight into the rate of rolling up of the 
trailing vortex sheet may be gained on the basis of some 
similarity considerations. Thus, consider two wings 
of similar span loading but differing in span 6 and total 
circulation Ip. The symbols referring to the refer- 
ence wing are denoted by asterisks; those of the second 
wing are plain. The ratio of the induced velocities g 
at corresponding stations behind the wings, neglecting 
the small influence of the difference in chordwise load- 
ings, is given by 


q/q* = (To/b)/(To*/b*) (3) 


The ratio of the distances d, in terms of wing spans 
from the trailing edge to stations having similar degrees 
of rolling up of the trailing vortex sheets is 

Vo q 


d b - = VobI'o* A CL 


= = = 4 
d*/b* * q* Vo*b*Ty A* Cc.” a 


where A and C, represent, respectively, aspect ratio and 
lift coefficient. From this condition one can conclude 
that the distance required for the trailing vortex sheet 
to become essentially rolled up is of the form 
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e/c = K(A/C,) (b/c) (5) 


where c represents the wing chord and A is, as yet, an 
unspecified constant. The point where the vortices 
are rolled up is hard to specify precisely, since the vor- 
tices only asymptotically approach the completely 
rolled-up condition as the distance from the wing ap. 
proaches infinity. To indicate the order of magnitude 
of this quantity, it will suffice to give the value A = 
0.28 for wings having elliptical span loading, found by 
Kaden"! in a lengthy theoretical study. This matter 
will be discussed later in this paper. 

We see from the foregoing equations that the degree 
of rolling up of the trailing vortex sheet depends not 
only on distance and lift coefficient but also on the geom- 
etry of the wing. As an illustrative example, compare 
the magnitudes of e/c for a rectangular wing of aspect 
ratio 6 and a triangular wing of aspect ratio 2, assum- 
ing both wings are twisted to produce elliptical load- 
ing. 


(e/c)rect, = 0.28 (36/C,), (e/c)tri. = 0.28 (2/Cz) (6) 


It is seen that the trailing vortex sheet rolls up 18 times 
more rapidly, in terms of chord lengths, behind the 
low-aspect-ratio triangular wing than behind the high- 
aspect-ratio rectangular wing. 
ratio, the vortex sheet may become essentially rolled 


For wings of low aspect 


up into two trailing vortices within a chord length or 
less of the trailing edge. 


Behavior at Stations Far Behind the Wing 


The behavior (i.e., the strength, position, motion, 
and core diameter) of the rolled-up vortices at stations 
far behind the wing is independent of Mach Number 
and can be determined by use of Eq. (2). First, how- 
ever, the strength and positions of the rolled-up vor- 
tices must be determined from three-dimensional con- 
siderations, since they are necessary to determine the 
boundary conditions. The strength, at either subsonic 
or supersonic speeds, of one of the trailing rolled-up 
vortices is equal to the sum of the strengths of all the 
vortices shed from one-half of the wing, and, hence, it 
is equal to the magnitude of the circulation I) around 
the wing in the plane of symmetry; thus, 


Myo = (Vo 2)(CiC) y =o (7) 


where c, represents the section lift coefficient. The dis- 
tance between the two rolled-up vortices 2s’ is specified 
by the fact that the lift impulse must be preserved 


throughout the rolling-up process; thus, 


25’ = ie (po Volo) = (C,S) (iC) y 0 (S) 
where Z and S represent, respectively, the lift and area 
of the wing and po is the free-stream density. The 
distance 2s’ has a simple geometric interpretation: Ona 


span-loading plot of T or qc versus y, the distance 2s’ 
represents the width of a rectangle having the same 


height and area as the span-loading curve. 
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TRAILING 


The further study of the asymptotic behavior at sta- 
tions far behind the wing may be made entirely on the 
basis of Eq. (2) taken together with the conditions ex- 
pressed by Eqs. (7) and (8). From these considera- 
tions, it can be determined that the two vortices move 
downward with a velocity 


saath Vo t7) Met )y a 0}? (C,S) (9) 


@. = —Vo/(4ms') = 


The actual displacements of the vortices from the xy 
plane cannot be ascertained at this point, however, 
since no statement can be made, as yet, of the downward 
This 


For 


velocity of the vortices at intermediate stations. 
question will be considered in the next section. 
elliptic span loading, Eqs. (7), (8), and (9) reduce to 
the following: 


C 4 G 
= “4 : (10) 


The rolled-up vortices far behind the wing are not 
idealized line vortices of potential theory, of course, 
but must have cores of finite diameter. A good working 
approximation at either subsonic or supersonic speeds 
may be had by simple energy considerations, assuming 
that the vortex cores are circular in shape and rotating 

The velocity potential of the flow 
outside of the vortex cores is 


I"y Z Z 
tan”! - — tan : (11) 
2a yo s ys 


The stream-line pattern corresponding to a vortex pair 
+s’ consists, as shown in Fig. 2, of two 


as solid bodies. 


situated at y = 
symmetrically arranged families of circles having radii 
+a’. The kinetic energy of 
the fluid outside of the vortex cores (per unit length in 


a and centers at +V 5” 

the stream direction) is given by 
, Po ¥ O¢ 

Eo = g tal dX — 

2 - on 


‘tvs? +a? — a, 
+ a.’ + a, 


log (12) 


2r ee 
where the integration is carried round the contour 
illustrated in the figure and a, represents the core radius. 
The kinetic energy of the fluid inside the vortex cores is 


simply 


. Po “a Ip fe a es poly” 
E, = 2 Qer dr = 
. ( ’) f Ea (5) ee Sir 


By equating the total kinetic energy per unit length to 
the induced or vortex drag D; of the wing, an expression 
may be found relating the radius of the core to the in- 


(13) 


duced drag, 


ol'o” 
Saw SR ow = 
St 
yf he We at ~~ 
1 + 4 log : (14) 


fs —VsP? +a? + a, 
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Since a, is much less than s’, Eq. (14) may be simplified 
to the form 


Po ‘ 


(1 +41 = 4) 15) 
og { re) 
St S a, 


The ratio of the core radius a, to the semispan of the 


D, = Cp; gs = 


vortices s’ is then given in terms of nondimensional 


parameters as 


a, twACp, (s'\2 177 
~=2 1 + exp —— ) - l (16) 
s Ci? \s 1) 
If the span loading is elliptical, Cp, = C,*/7A and 
s'/s = w/4and the core radius is given by 
a. = 0.197%s’ = 0.155s (17) 


The Rolling Up of the Vortex Sheet 


The complete three-dimensional determination of 
the shape of the trailing vortex sheet throughout the 
rolling-up process presents a problem of extreme diffi- 
culty. Examination of the velocity field in the wake 
on the basis of Eq. (5), however, reveals that in most 
practical problems it is possible to introduce the im- 


portant simplifying approximation of treating the 
motion of the trailing vortices by means of Eq. (2) 
rather than Eq. (1). Generally speaking, Eq. (2) 


yields good results when applied to flow fields that 
change slowly in the stream direction and is directly 
applicable to the present problem at all stations behind 
low-aspect-ratio pointed wings. This is the most im- 
portant case, since, as indicated by Eq. (5), it is only 
with low-aspect-ratio wings that the trailing vortex 
sheet is likely to be essentially rolled up at the stations 
customarily occupied by the tail. Although Eq. (2) 
obviously cannot be used to determine the forces on 
high-aspect-ratio straight wings, it might still be ex- 
pected to yield useful results regarding the relative mo- 
tion of the vortices behind such wings, since Eq. (5) 
indicates in these cases that the rolling-up process pro- 
ceeds at a slow rate. On the basis of the foregoing 
statements, the problem of determining the shape of the 
trailing vortex sheet in spatial compressible flow be- 
comes analogous to the problem of determining the 
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motion of a distribution of free vortices in a two- 
dimensional nonsteady incompressible flow. 

Considerable simplification may be achieved in the 
present problem by considering the motion of the center 
of gravity of the vortex system rather than by tracing 
the motion of the individual vortices. While generally 
useful, such considerations are particularly valuable in 
the establishment of the position of the fully rolled-up 
vortices far behind the wing. The motion of the center 
of gravity of a vortex system may be determined by 
several means; probably the simplest procedure for 
the present application is that of Betz.'” 

Consider a group of vortices in a two-dimensional 
region bounded by rigid walls. If the 
held rigidly in place, there exists on each vortex a force 


vortices are 


having components R;, = —pow,'; and R;, = pov ly; 
as given by the Kutta-Joukowsky law, where I’; repre- 
sents the circulation of the vortex, clockwise rotation 
being considered positive, and where v; and w, represent 
the lateral and vertical velocity components at the 
position of the vortex were the vortex not there. The 
resultant of the forces on all the vortices must be equal 
to the negative of the resultant force exerted on the 
boundary walls; thus, R, = —2R;, and R, = —2R;.. 
If the vortices are now set free, the center of gravity of 


ly 


the vortex system moves, relative to the rigid walls, 
with velocity components 7 and @ given by the rela- 


tions po®2T, = porw l; = R, and pwrl; = —R, or by 
w= R, pox’ ;, v= —R, pol’, (18) 


The foregoing theory will be applied to the present 
problem by introducing a rigid wall enclosing one-half 
of the vortex system as shown in Fig. 3. If the en- 
closed region is allowed to grow in size, the contribution 
to the resultant force given by the pressures on Segment 

+ The center of gravity of a group of vortices is defined in a 
manner analogous to the definition of the center of gravity of a 
The coordinates 7 and 2 of the 


zr; 


similar field of point masses. 
center of gravity of the vortex system are given by 7 = Ly,I'; 


and § = Ys,I':/ZTs. 
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I] of the boundary approaches zero as the boundary 
approaches infinity. Since the resultant force is always 
directed perpendicular to Segment I lying in the plane 
of symmetry, R, must equal zero. Consequently, the 
lateral position of the center of gravity of the vortices 
trailing from each half of the wing remains invariant 
throughout the rolling-up process and may be deter- 
mined most readily by considering the c.g. position of 


the flat vortex sheet 


” Syl; l s aT 
Pf 6 See ae en y dy = 
wij; Ip 0 dy 


The vertical velocity @ of the center of gravity of the 
vortices trailing from each half of the wing is found from 
Eq. (18) after substituting the proper expressions for 
R,and =V;. For the right half wing, R, is given by 


i Po -— 
= [ (p — po) dz = + : f (wy —0)* dz 


The center of gravity of the vor- 


const. (19) 


and =I; equals —T». 
tices at any station x therefore moves vertically with a 


| oo oe: 
W(x) = — — / (w(x), ~0]* dz 


av 


velocity 
(20) 
Immediately behind the trailing edge, where the vor- 


tices lie in a straight line as shown in Fig. 4, the induced 
velocity at the wall is given by 


# | 
T JO 


The velocity of the center of gravity of the free vor- 


| Pa - 
aa / dz | / 
1*1"p 0 J 0 


As a more specific example, consider the case of 
If the vortices are held fixed, the 


yi(dl/dy,)dyy 
(yy? "+ 3°) 


Ww (Xp E i e = 


tices is then 


W(X E ) =—__ (21 ) 


yi(d | 
(yi? + 27) 


elliptic span loading. 
flow is similar to that around a flat plate moving down- 
ward with the velocity 








w= —Iy/b = —2V0C,/rA 
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TRAILING VORTEX SHEE 


If the vortices are set free, the velocity of the center of 
gravity of the vortices decreases suddenly to 


W(X E.) = —(1 ee 4) (To b) = — 0.43 VoC, WA (22) 
Far behind the wing, where the vortex sheet is con- 
sidered to be fully rolled up into two line vortices, the 
corresponding velocity is 


w( oo ) = —27)/r*b = —0.41 VoCL TA (23) 
It is seen that the center of gravity of the trailing vortex 
sheet moves with nearly the same velocity at the two 
limiting stations. 

A more detailed analysis of the behavior of vortex 
systems can always be made on the basis of Eq. (2) 
by replacing the continuous sheet of vortices with a 
finite number of discrete line vortices and calculating 
step by step the motions of each vortex. Such a cal- 
culation has already .been carried out for the case of 
elliptic loading.'* In this particular treatment, the con- 
tinuous vortex sheet was replaced by 20 vortices of 
equal strength, and the results were presented by giv- 
ing, both numerically and graphically, the positions 
of each of the vortices at different times. 
These results have been applied to the wing problem 
by relating the time to the distance behind the wing. 
Fig. 5 shows the positions of the vortices at several sta- 


several 


tions behind the wing. 

It may be seen that the vortex cores extend down- 
stream in nearly the free-stream direction, whereas 
the center of the vortex sheet becomes displaced sub- 
stantial distances below the xy plane. The variation 
with distance of the positions of the center of the vortex 
sheet and of the vortex center of gravity is shown in 
Fig. 6. As is indicated in this figure, the position of 
the vortex center of gravity may be prescribed within 
narrow limits by assuming that it moves with a con- 
stant velocity given first by Eq. (22) for the flat vortex 
sheet and then by Eq. (23) for the fully rolled-up vor- 
tices. Also shown in Fig. 6 are the limiting positions 
for the center of the vortex sheet as calculated in a 
similar manner. 

Since the discrete vortex method just 
lacks generality because of the necessarily nonanalytic 


described 


character of the solution, results of an analytical in- 
vestigation by Kaden!'! of the behavior of a continuous 
vortex sheet, according to Eq. (2), will be reviewed and 
applied to the present problem. The precise problem 
solved in reference 11 is that of the rolling up of a vor- 
tex sheet of semi-infinite width, where the initial vor- 
ticity corresponds to steady flow around the edge of a 
two-dimensional flat plate of semi-infinite lateral ex- 
tent. Although the results are strictly applicable to a 
vortex sheet of finite width during only the early stages 
of the rolling-up process, they represent a solution that 
is useful in the formulation of some further approxima 
tions to the paths of the vortex cores and to the distance 
required for the trailing vortices to become essentially 


rolled up. 
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Fic. 5. Shape of vortex sheet for elliptic loading 
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The results of this investigation may be expressed 
most conveniently by introducing a new coordinate 
system (Fig. 7) such that the 7 axis coincides with the 
position of the flat vortex sheet and the origin lies at 
the edge of the sheet. In applying these results to the 
wing problem, it is important to note that this coordi- 
nate system is not fixed with respect to the y and 2 
axes but moves downward with the center of the vortex 


sheet. If the circulation distribution in the vicinity of 
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the edge of the flat vortex sheet is assumed to be of the 
form I! = 20~/n, the coordinates of the centers of the 


vortex cores are 


ne = 0.340(0d/Vo), £- = 0.525 (od/ Vo)? (24) 


It will be noted that y,, as given by the above equation, 
does not have an asymptote corresponding to the lateral 
position of the fully rolled-up vortex. This is a conse- 
quence of replacing the finite span with a semi-infinite 
one. The trailing vortices are considered to be essen- 
tially rolled up at the distance e behind the wing, where 


ne is equal to (s — s’) as indicated in Fig. 8. 





a 

















Fic. 8. 


Eq. (24) may be generalized slightly and put in non- 
dimensional form by letting 


o = (1/2) lim P/V = (Vo/4) lim (cc/V n) 


7 —~ 0 n— 0 


The coordinates of the centers of the vortex cores are 
then given by 


ne = (s — 8’) (d/e)”, f= 1.54 (s — s’)(d/e)” (25) 
and the distance e by 
e = 20.2 (s — s’)’*/lim (qc/V n) (26) 


nO 
For elliptic span loading, the foregoing equations re- 


duce to simple expressions. Thus, since the span load- 


ing is represented by 
ac = (2C,S/mrs) V1 — (y/s)? = 
(2C,S/ms) V 2(n/s) — (n/s)? 


the distance e to the point where the trailing vortices 
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are considered to be essentially rolled up is expressed by 
e/c = 0.28(4/C,)(b/c) (27) 


The coordinates of the centers of the vortex cores are 
then 


_ ee 5 ¢ C, d 
> os oe “) tata ( ? 3 (28) 
Ss  . 4 S Ss Z 4 Ss 


Eq. (27) should be compared with Eq. (5) derived by 
similarity considerations. If the span loading is other 
than elliptical, the rolling up, of course, will proceed at 
a rate other than that just indicated. In general, 
the vortices will roll up faster if more of the load is con- 
centrated near the wing tips (as occurs with rectangular 
wings at subsonic speeds) and conversely. To illustrate 
the importance of this effect, calculations have been 
made for the distance e/c for rectangular wings of vari- 
ous aspect ratios. These calculations were based on 
theoretical span load distributions calculated by Glau- 
ert® by means of Prandtl’s lifting-line theory. The 
distance eis then given by 


ajo = E'/C, (29) 


where A’ is a constant for any given wing depending 
upon the plan form and the span load distribution. 
The values of AK’ are plotted as a function of aspect 
ratio in Fig. 9, together with the corresponding values 
for rectangular, elliptical, and triangular wings twisted 
so as to produce elliptical span loading. In every case, 
the chord c implied is the root chord. 

As previously pointed out, Kaden’s theory represents 
a good approximation for the lateral position of the 
vortex core at stations close behind the wing but fails 
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a/e = .09 a/c = .35 





Photographs of wake at various stations behind a 


Fic. 10 
a = 12°, CL = 0.55. 


triangular wing of aspect ratio 2. 


to approach the known asymptotic position far behind 
the wing. To bridge this gap, an empirical interpola 


tion formula is proposed 
ne = (s — s’) tanh (d/e)” (30) 
which coincides at small distances with Kaden’s for- 
mula and, at great distances, approaches the known 
asymptote. 
Experimental Studies.— The 
vortex systems was studied further by means of simple 
These experi- 


behavior of | trailing 
visual-flow experiments in a water tank. 
ments were conducted by driving a model wing ver- 
tically into the tank at a constant velocity and photo- 
graphing the water surface with a movie camera. The 
trace of the trailing vortex sheet was made visible by 
applying fine aluminum powder to the trailing edge of 
the wing before each run. The distance of the model 
below the surface was indicated on each photograph 
by the position of a moving tape. The models were of 
8-in. span and had flat-plate profiles, except for rounded 
leading edges and sharpened trailing edges. Although 
a large number of wings having various aspect ratios and 
plan forms were tested, only the results for a triangu 
lar wing having an aspect ratio of 2 and an elliptical 
wing having an aspect ratio of 5 will be presented here. 

The results of each run were recorded by a series of 
approximately 40 photographs. Abridged series of 
photographs are shown in Figs. 10 and 11 for the tri- 
angular wing at angles of attack of 12° and 20°. The 
projections in the free-stream direction of the wing-tip 
positions are indicated in the photographs by the inter- 


ife = 1.80 


a/e = 1.45 















Fic. 11. Photographs of wake at various stations behind a tri- 
angular wing of aspect ratio 2. a = 20°, Cr 0.93. 
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Fic. 12. Water-tank measurements of the lateral position of the 


vortex cores behind an elliptical wing of aspect ratio 5 


sections of the vertical lines and the horizontal 
markers. Although it is not the purpose to obtain 
exact quantitative information from these experiments, 
the photographs clearly illustrate the rapidity with 
which the trailing vortex sheet rolls up behind a low- 
aspect-ratio wing and verify, at least qualitatively, the 
predictions of the theory regarding the shape and loca 
tion of the vortex sheet, its rate of rolling up, and the size 
of the vortex cores. 

To illustrate further these results, Figs. 12 and 15 
have been prepared showing the paths in plan view of 
the vortex cores behind the elliptical and triangular wing. 
The results for the elliptical wing show that the vortex 
cores approach the asymptotic spacing for elliptical 
loading (s’/s = 0.785) at stations several chord lengths 
behind the wing. With the triangular wing, the vor 
tex cores approach their asymptotic spacing more 
rapidly, as is indicated by the theory, but approach a 
different asymptote at each angle of attack. The 
variation of the asymptotic spacing cannot be pre- 
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Fic. 13. Water-tank measurements of the lateral position of the 


vortex cores behind a triangular wing of aspect ratio 2. 
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Wind-tunnel measurements of the span loading on a 
triangular wing of aspect ratio 2.04. 
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Fic. 15. Asymptotic spacing of rolled-up vortices far behind 


triangular wing of aspect ratio 2. 


dicted on the basis of the theoretical span loading" but 
can be explained on the basis of the variation of the 
experimental span load distribution. Fig. 14 shows 
plots of the span load distribution measured in a wind 
tunnel on a 25-ft.-span model of a thin triangular 
wing having an aspect ratio of 2.04 and a modified 
double-wedge profile with a rounded leading edge. 
The asymptotic spacing of the vortex cores was cal- 
culated from each of the span loading diagrams by 
means of Eq. (8) and was plotted as a function of the 
angle of attack in Fig. 15, together with the values 
indicated in Fig. 13. The curve for the experimental 
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location is shown dashed between 8° and 12° because 
the vortices were still moving together at the last sta- 
tion behind the wing. It is seen that these consider 
ations adequately account for the observed behavior. 


DOWNWASH AROUND A SWEPT HORSESHOE VORTEX 


The fundamental role of the horeshoe vortex in classi- 
cal wing theory is well known. 
of this familiar concept to the general case, the down- 


To enable the extension 


wash field of a swept horseshoe vortex of arbitrary 
sweep lying in the xy plane (see Fig. 16) has been deter- 
mined for both subsonic and supersonic speeds. In this 
task it is necessary, because of the essentially three- 
dimensional character of the problem, to work with 
Eq. (1). 

To utilize existing concepts of supersonic, as well as 
subsonic, flow theory, the induced velocity field of the 
swept horseshoe vortex will be determined by consider- 
ing the flow about an equivalent doublet sheet." 
Such a doublet sheet represents a surface of discon- 
tinuity in the xy plane across which there is a jump 
Ags in the value of the velocity potential equal to the 
magnitude of the circulation I of the vortex. The 
lift of this vortex system is given at both subsonic and 
supersonic speeds by the Kutta-Joukowsky law 


2Vopo J 0 Ags = 2Yopo I ol 


where yo is the semispan of the horseshoe vortex. 

The general solutions of such boundary-value prob- 
lems were given in references 15 and 16. For subsonic 
flow the solution is 


1 Oo a ‘oO 1 
/ / Ags ( dx, dy, (31) 
iwiOs J, . 02; p/s ; 


W(X, y, 2) = 


where 











Swept horseshoe vortex 
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and 8 = V1 — M*. The integration is carried over and B = VM? — 1. The integration for the super- 
the portion 7 of the xy plane enclosed within the modi- sonic case is carried only over the portion 7 of the above 
fied horseshoe vortex. For supersonic flow, the solu- area which lies forward of the trace in the xy plane of 

the Mach forecone with vertex at the point (x, y, 2). 
— The sign [~ is to be read “the finite part of’’ and has 


1a se a 1 
w(x, ¥, 2) = Ag: ( \ dx, dy, (32) _ the property 
we 2m dz tf a OS 


[ f(x) j [ f(x) — f(b) j 2f(b) 
shere dx = ax — 
where Je (6 — x) Je 6-2) Vb-a 


a 


p = V(x - x)? _— B2[(v - \1)? + + ie 2)? | (3. 


tion 1s 


When the operations indicated in Eq. (31) are performed, the downwash corresponding to a subsonic swept 


horseshoe vortex of circulation I is given for any point (x, y, 2) in space by 


sr +> | x — kBmy }} mx + kBy 4 
= far (x — kBmy)? + B?27(1 + m*)J la/ x2 + B%(y? + 22) 


k=-l , 
(1 + m*)Byvo — (mx + RBy) ( 4 =. | ky + Yo | Y 
fe a Bm)? +- B?[(y = ky)? + 2]! alk (y oe kyo)? + 2? 


k +1 
—1 


k= +1 
{ x — Bmyo \ : 
)! + (34) 
V (x — Bmyo)* + B?[(yo + Ry)? + 2°]! 
In applying Eq. (34), it is important to remember that infinite values of downwash are predicted along the entire 


vortex line. 
In a similar manner, Eq. (32) gives the downwash at any point in space due to a supersonic swept horseshoe vor- 


tex. 
ér | x— kB my | j mx — kBy 
w=r.p./ — = = 
2m > (x — kBmy)? — B22(1 — m?) SIV x2 — §2(y? + 22) 
fo + 
(m? — 1)Byo — (mx — kBy) \ r 7 a (yo + ky) (x — B myo) 35) 
= = ae = = (39 
V (x — Bmyo)? — B*[(yo — ky)? + 22] ” — [(yo + Ry)? + 22]V (x — Bmyo)* — B?|(vo + ky)? + 27] 
k +1 


where r.p. signifies the real part of the bracketed quantity. For supersonic flow, Eq. (35) indicates infinite values 
for the downwash at all points on the Mach cones extending downstream from the three corners of the swept 
horseshoe vortex (except in the xy plane), as well as along the vortex line itself. 


DOWNWASH AROUND WINGS The bound vortices should, in general, be concentrated 


. a ; — in a lifting line lying as close as possible to the locus 
It is a well-known principle in subsonic wing theory 8 rs —— the locu 


that the induced velocities produced by a thin lifting 
wing may be calculated by considering a chordwise 
and spanwise distribution of vortices. The flow field 
of a wing at supersonic speeds can be built up in an 
manner."’~'® In general, these vortices position of the lifting line can always be determined if 


of the chordwise center-of-pressure positions. Thus, at 
subsonic speeds, the lifting line is generally located 
along the locus of the quarter-chord points. At super- 
sonic speeds, no such statement can be made, but the 


analogous 
would be distributed chordwise and spanwise in such a _ the load distribution is given. 
manner that the specified load distribution would be 
satisfied at every point of the wing plan form. Such 
calculations may be simplified, however, since, as indi- 
cated in references 10, 16, and 19, the error in the down- 
wash introduced by neglecting the chordwise distribu- 
tion of the bound vortices is negligible at points beyond 
a fraction of a chord length from the trailing edge. generally, by Eq. (26). 


The following discussion of the downwash behind 
wings will be divided into sections depending upon the 
degree of rolling up of the vortices as indicated by the 
relative magnitudes of the distance from the trailing 
edge d and the distance e given by Eq. (5) or, more 
The important special cases 
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where d << e and where d >> e will be considered 
first; finally, the intermediate case will be discussed. 


DOWNWASH AT d << e 


If dis much less than e, the downwash behind a wing 
traveling at either subsonic or supersonic speeds may 
be approximated by assuming that the trailing vortices 
lie in a flat sheet extending downstream from the trail- 
ing edge. This assumption is usually satisfactory for 
the calculation of the downwash at the tail of conven- 
tional subsonic airplanes and forms the basis of most 
downwash studies. In general, the downwash velocity 
at any point may be determined by superposing a dis- 
tribution of swept horseshoe vortices in accordance with 
the given span loading. In this manner, an expression 
for the downwash velocity in subsonic or supersonic 
flow may be determined from Eqs. (34) or (35), respec- 
tively, by considering yo as a variable y,, replacing the 
circulation [ with —(dI'‘/dy,)dy,, and integrating the 
resulting expression from y; = —s toy, = +s. This 
integration, however, is extremely complex, in general, 
although many important special cases for both sub- 
sonic and supersonic speeds have been presented in the 
literature. 

It is often more practical to calculate the downwash 
velocity by superposing a finite number of horseshoe 
vortices sswpt in accord with the geometry of the wing 
and distributed so that the span loading is approxi- 
mated in a stepwise manner. Such a treatment, which 
was presented in great detail for unswept wings in in- 
compressible flow in references 10 and 20, proceeds in 
the following manner. The span load distribution is 
approximated by an equivalent stepwise distribution, 


keeping the area under the two curves equal. A vortex 
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of strength given by the amount of the rise is considered 
to be shed at each step. The downwash velocities asso- 
ciated with each component vortex are calculated by 
means of Eqs. (34) or (35) and summed to obtain the 
total value. Care must be taken to avoid attempting 
to evaluate the downwash at points on, or very near, 
either the trailing vortices or, except in the xy plane, 
on the Mach cones originating at the corners of the 
horseshoe vortices. On these lines and surfaces, in- 
finite velocities are introduced purely as a result of the 
replacement of a continuous distribution of vortices 
with a stepwise approximation. 

The calculation of the downwash field 
greatly simplified in the regions where Eq. (2) is a good 
approximation to Eq. (1). In the cases where the lift 
coefficient is sufficiently small that d < < e, the down- 
wash at any Mach Number is 


I [- sdY (yw, — y)dyy 
21 _ (36) 
2r J-; dy [(m — y)? + 27] 


For elliptic loading, the solution of Eq. (36) is given in 


becomes 


wy, z) = — 


many places.®?! 


sinh uw cosh pu oad 
(37) 


2Cz 
w - mci 9 . 9 
TA cosh? w — sin? X 


where y = } cosh uw cos \ and zg = bsinhywsind. A 
plot of these results is shown in Fig. 17. 

The accuracy of the theoretical downwash results 
may be improved by considering the trailing vortices 
as being displaced from the xy plane. To a first ap- 
proximation, the differential equation of the path of a 
vortex filament trailing behind a wing is dx = (Vo/w- =o) 
dz, which, when integrated, gives the vertical displace- 
ment of each vortex filament 


as = / (W:—0 Vo)dx (38) 
J ©T RK, 


If it is desired to calculate the downwash only in the 
plane of symmetry, the effects of the lateral distortion 
of the vortex sheet are small, and the entire vortex 
sheet may be considered to be displaced by the amount 
calculated at the centerline. 

For regions outside of the plane of symmetry, it is 
apparent that this last assumption cannot be particu- 
larly accurate and that recourse must be had to more 
accurate methods. One method of doing this would 
be to determine the shape of the vortex sheet by means 
of Eq. (38) and to calculate the downwash anew with 
the vortices displaced into their new positions. This 
method is not entirely satisfactory, however, when 
working with a continuous distribution of vortices, 
since the displacement of the vortex sheet from the xy 
plane greatly complicates an already difficult integra- 
tion. An alternative procedure based on distributing 
double vortex sheets in the xy plane to simulate the 
sinking of the original vortex sheet has been proposed by 


von Karman and Burgers,*! but the mathematical 
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TRAILING VORTEX 
complications are comparable with those encountered 
with the distorted vortex sheet. 

It must be noted that the procedures just described 
for obtaining the correction for the displacement of the 
trailing vortices can only be used in cases where the 
distortion of the vortex sheet is small. These methods, 
therefore, are of no avail for the investigation of the 
downwash in regions where the vortices have become 
rolled up to any appreciable degree. This problem 
will be discussed in a later section. 


DOWNWASH AT d >> e 


At stations where d is much larger than e, the in 
duced flow field of a wing may be satisfactorily ap 
proximated by replacing the wing with a single swept 
horseshoe vortex placed as shown in Fig. 18. Accord 
ingly, the downwash may be computed by using Eqs. 
(34) or (35) for the downwash field of a swept horseshoe 
vortex lying in the xy plane and applying corrections 
for the vertical displacement of the vortex sheet and for 
the finite size of the vortex cores. The strength and 
span of the horseshoe vortex are given by Eqs. (7) and 

8), respectively, while the location and sweep of the 
bound vortex are determined as described previously at 
the beginning of the discussion of the downwash around 
wings. 

The correction for the sinking of the trailing vortices 
may be approximated by inclining the horseshoe vortex 
system away from the xy plane so that each of the trail 
ing vortices coincides as closely as possible with the 
center-of-gravity line of the vortices trailing from the 
corresponding half of the wing. Since it was shown that 
the center-of-gravity line is inclined at nearly the same 
angle at all stations behind the wing, a suitable ap 
proximation may be made in most cases by using Eq. (9) 
for the downward velocity of the vortices at great dis 
tances behind the wing. Thus, the plane of the equiva 
lent vortex system is inclined at angle 6 with the xy 
plane such that 


(59) 


6 = 0 imC.S 


— (Ci ), 


To a first approximation, the downwash may be calcu 
lated by using Eqs. (34) or (35) for a swept horseshoe 
vortex lying in the xy plane and by replacing z in the 
equation with z’ = s — 6 (x — xp). The downwash 
within the vortex cores may be approximated by as 
suming that each core has a radius given by Eq. (17) 
and rotates as a solid body. 

In most actual applications where d >> e, Eq. (2) 
represents a good approximation to Eq. (1). In these 
cases the downwash field can be computed at once for 
all wings. The results of these calculations are shown 
in Fig. 19 as a function of y/s’ and the vertical elevation 
A comparison 


, 


z’|/s’ from the plane of the vortices. 
of these results with those of Fig. 17 illustrates the 
marked effect on the downwash produced by the rolling 
up of the trailing vortex sheet. 
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Equivalent vortex system for downwash calculation at 
stations where d > > e. 
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Fic. 19. Downwash far behind wing for completely rolled-up 
vortices 


DOWNWASH AT d = e 


The remaining problem of calculating the down 
wash for regions where d is of the same order 
of magnitude as e is the most difficult problem of 
all. The difficulty lies not only in the determination 
of the vortex paths but also in the integration of the 
resulting expressions for the downwash. Fortunately, 
a large number of practical problems fall into one or the 
other of the two preceding cases so that oftentimes 
one is not concerned with this phase of the downwash 


problem. 


If the vortex sheet rolls up slowly behind the wing 
so that Eq. (2) may be used instead of Eq. (1), the 
downwash may be calculated by determining the posi- 
tion of a finite number of vortices in a step-by-step 
process and then calculating the associated downwash 
velocities. If the span loading is elliptical, the vortex 
positions shown in Fig. 5'* may be used directly. 
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Fic. 20. Various approximations to the lateral position of the 


vortex cores. 


If the conditions are such that Eq. (1) must be used, 
the problem becomes extremely difficult. The only 
attempt at a solution for a problem of this type is that 
of Lotz and Fabricius,”? in which the downwash on the 
x axis produced by straight wings in incompressible 
flow is determined. Since this method involves several 
novel elements and can be extended to other applica- 
tions, their approach to the problem will be reviewed. 


It is first assumed that the vortices all lie in the xy 
plarie so that the rolling-up process is considered only 
through the lateral displacements of the vortices. 
Such an approximation is usually satisfactory for re- 
gions close to the plane of symmetry but would have 
to be abandoned in the general vicinity of the vortex 
Further, the paths of the partially rolled-up 
vortex cores are assumed, on the basis of an arbitrary 
modification of the theoretical paths of Kaden, to 
The shape of this 


cores. 


be given by 9, = (s — s’) tanh (d/e). 
curve and its significance as an approximation on the 
basis of a straight-line replacement of Kaden’s vortex 
system are clear from Fig. 20. Since Kaden’s theory is 
most accurate at stations near the wing, this approxi- 
mation could be improved upon without unduly com- 
plicating the calculations by using a curve such as given 
by Eq. (30) which is tangent to Kaden’s curve near 
the wing. Third, it is assumed that the strength of 
the vortices lying between the cores is given by 


(or, Oy), = (Or OY) + g.° P(x) 


where 0 < P< 1. With the foregoing assumptions, 
the strength I’, of the vortex cores at any station is de- 
termined, since the total strength I) of the vortices 
trailing from each half of the wing remains constant 
at all stations; thus 


r. = Ty — P(x) [fo — Py.) ] 


Further, P(x) is fixed by the requirement that the 
lateral position of the center of gravity of the vortices 
trailing from one side of the wing is invariant 


y(x) — fe ly )dy 
0 


P(x) = - — 


*ye(x) 
V(X) = I (T I) )dy 
J 0 


y(x) —s 


ye(x) 
v(x) — / (T'/To)dy 
J 0 


(40) 


SCIENCES—JANUARY, 1951 

One sees that P(x) takes the value P(x; ,) = | at the 
trailing edge of the wing, since here y,(x) = s; whereas 
at great distances behind the wing P(~) = 0, since 
y-(x) approaches s’ as x approaches infinity. 

By means of the foregoing assumptions and equa- 
tions, the vorticity distribution is determined through. 
out the entire vortex wake. The downwash may then 
be calculated for subsonic or supersonic speeds by using 
Eqs. (31) or (32) for the downwash due to an equivalent 
doublet sheet. The doublet intensity is specified on 
the wing as before and in the wake by means of the fore- 
going statements. For points along the x axis, Eqs 
(31) and (32) simplify, respectively, to the following 


form: 


] iy . Agsdx,dy, 
w(x, 0,0) = / 5 a | 
4a Je. (x — %)* + 6,7)” 


0. 0) l ae Agsdxidy, 19 
W(x, UV, = = 4 (+42 
27 ; / [(x — %)* sins By," | : 


where the doublet intensity is given by 
Ags = l'r.2.(m)°P(%1) + To — P(m) [To — T(y-)] (48 


Despite the assumptions and simplifications of the 
foregoing analysis, the calculations still prove to be 
extremely difficult. In the work of Lotz and Fabricius, 
the theory was further restricted to unswept wings in 
incompressible flow; even then it was necessary to re- 
sort to graphical or numerical methods of integration 
Attempts to extend this theory to regions away from 
the x axis or to eliminate some of the restrictions have 
been discouraged by the difficulty of the integrations. 
The labor required to carry out the calculations simply 
appears to be incommensurable with the value of the 
result, considering the doubtfulness of some of the 
simplifying assumptions. This conclusion is further 
amplified when one considers the problem of predicting 
the downwash at the tail of an airplane where the ideal- 
ization involved in considering the wing alone is already 
rather extreme. 
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Gust Loads on Rigid Airplanes with Pitching 
Neglected’ 


R. L. BISPLINGHOFF,?t G. ISAKSON, fanp T. F. O’BRIEN** 
Massachusetts Institute of Technology 


SUMMARY 


The motion of a rigid airplane flying into a linearly graded gust 
The pitching motion of the airplane is neglected, 
The 
The 


equation of motion is solved by application of the Laplace trans- 


is analyzed. 
but unsteady aerodynamic effects are taken into account. 
analysis is first carried out assuming incompressible flow 
formation. The solution is expressed in the form of an allevi- 
ation factor and is found to depend upon two nondimensional 
parameters, one being the so-called mass parameter and the other 
Com- 
prehensive data on the alleviation factor are presented for a 
practical range of values of the mass parameter and the gradient 
The case of flight at high subsonic speeds is considered, 


relating the gust gradient distance to the airplane chord. 


distance. 
and it is found that an analysis of this case is not immediately 
practical, since a great deal of basic numerical data on the aero- 
dynamic forces is required and little of these data has yet been 
computed. On the other hand, it is found that the case of super- 
sonic flight may be readily analyzed, and this is done, under 
the assumption of two-dimensional flow, for flight through a 
sharp-edged gust. Alleviation factor data are presented for a 
number of values of the Mach Number and for a practical range 
of values of the relative density parameter. It is shown that gust 
alleviation due to unsteady flow is likely to be slight for any 


practical supersonic airplane 


List oF COMMON SYMBOLS 


p = air density 

l = airplane horizontal velocity 

S = wing area 

dC, /da = slope of lift curve for airplane 

Ww = gust velocity 

m = mass of airplane 

= vertical displacement of airplane center of gravity 

positive down 

p = Laplacian operator 

r = airplane dimensionless mass parameter, 
m p(C 2)S(dCL, da) 

t = time 

s = dimensionless variable, (2U/c)t 

r = wing mean aerodynamic chord 

M = Mach Number 

iv = airplane relative density, m/p(c/2)S 


INTRODUCTION 
, | Nie SPECIFICATION OF GUST LOADS in current air- 
worthiness requirements is based upon an empirical 
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approach that has remained unchanged for a good many 
years. This approach consists of assuming that the 
effect of the gust is simply that of a sudden increase in 
lift on the airplane equal to (1/2)pl*S(dC,/da)(w/ U) 
or the increase in lift that would be obtained under 
steady flow conditions as the result of a change in angle 
of attack of w/l’.' It is seen that this is the result that 
would be obtained if the airplane were to fly into an 
idealized sharp-edged gust of velocity w and if it could 
be assumed that the entire airplane were instantly en- 
gulfed by the gust and that aerodynamic forces were 
instantaneously adjusted to steady-state values. 

In reality, these circumstances are never encountered. 
A sharp-edged gust does not exist. The entire airplane 
is not engulfed instantly by the gust, since a finite 
amount of time is required for the gust to progress over 
the airplane, and aerodynamic forces do not adjust in- 
stantaneously to steady-state values when there is a 
rapid variation of the angle of attack. The velocity w 
must therefore be regarded as that of an equivalent 
sharp-edged gust that will produce the same normal 
acceleration of the airplane under idealized conditions 
as the actual gust produces under real conditions. 
Specified values for w have been determined on this 
data normal accelerations 


basis from statistical on 


measured in flight through gusty air.” 


The relation between the maximum velocity in the 
actual gust and the velocity of the equivalent gust de- 
pends to a large extent upon the characteristics of the 
airplane itself, since the maximum lift does not occur 
until the airplane has penetrated some distance into the 
gust and has attained velocity in the direction of the 
gust. It is therefore necessary to take some account of 
this effect in determining loads from the sharp-edged 
gust formula. This is done by including an “‘allevi- 
ation”’ factor that is considered to be a function of the 
wing loading. It is, in reality, a function of additional 
parameters such as wing mean chord, the slope of the 
airplane lift curve, and parameters determining the 
pitching motion of the airplane. 


Numerous investigations into the effect of these 
additional parameters have been carried out for the case 
of a rigid airplane in low-speed flight. The unsteady 
flow effect has been studied by Rhode,’ Kiissner,* and 
Greidanus‘ for the case when pitching is neglected. 
The pitching effect has been considered by Bryant and 
Jones,® assuming steady flow; by Greidanus;! and by 
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Response to a linearly graded gust using steady flow 
theory 


Fic. 2. 


Greidanus and van de Vooren,® taking unsteady flow 
into account. 

The present work is similar to that of Greidanus' for 
the “‘pitching neglected”’ case, except that a gust with 
linear gradient is used rather than the cosine form 
assumed by Greidanus. In addition, consideration is 
given to the problem of gust loads for flight in the high 
subsonic and supersonic speed ranges. 


Gust Loaps IN LOoW-SPEED FLIGHT 


The first case to be considered is that of flight at a 
speed that is sufficiently low so that the air may be 
assumed incompressible. The following assumptions 
are made: 

(a) The airplane is initially in horizontal flight at 
constant velocity U. 

(b) The gust is blowing vertically and has a velocity 
distribution w in the direction of flight and a uniform 
velocity distribution in the spanwise direction. 

(c) Variation in the forward velocity of the airplane 
as it traverses the gust can be neglected. 

(d) The entire lift on the airplane may be applied to 
the wing alone in determining the alleviation factor. 

With the effect of the pitching motion neglected, the 
airplane consequently has but a single degree of free- 
dom, and its disturbed motion is governed by the 


differential equation 


—L (1) 


mz = 


where m is the mass of the airplane; 2 is the vertical dis- 
placement of the airplane’s center of gravity from its 
initial position, positive down; 2 is the second derivative 
of z with respect to time; and L is the lift, positive up. 


JANUARY, 1951 


Assuming first that the flow is steady and the chord- 
wise distribution of gust velocity on the wing is constant 
at the value that it really has at the wing mid-chord, 
the lift is given by the expression 


l 2 d¢ L (° g ) 
L= U2S a (2 
gen ee We” . 


where p is the air density, S is the wing area, dC, /da is 
the slope of the lift coefficient curve for the airplane, and 
2 is the first derivative of z with respect to time. 


The equation of motion can now be written as 


l [ 26 d¢ L & : g ) 
—-pU. + 3 
f ~ dal’ U 


Applying the Laplace transformation to this equa- 


mz = 


tion and using the symbol p to denote the Laplacian 


operator, we obtain 
° l ‘ dC, 
mp*s(p) = — pUS ; 


ada 


[@(p) ~ pz(p)] | 


where 2(p) and @(p) represent the Laplace transforms of 
z and w, respectively. 
Solving for 2(p), 
Lay W@(p) 
a(p) = —- : (2) 
Ne plp + A/A)(U/o)] 
where c is the wing mean aerodynamic chord and X is a 
dimensionless mass parameter defined as follows: 


A = m/p(c/2)S(dC,/da) (6 


Applying Heaviside’s partial fractions expansion and 
the convolution integral in the inverse transformation, 
the following solutions are obtained: 


*t 
1/r) (Ue T _ 
sities / le ‘ _ l |z0( T) drt (4 
J 0 
and 
1 vu? f*' LU 
oe —(1/A) (U/ce) 0 T 
. == — é wir) dr — w(t) 
vc 0 Ac 


(S 


It is seefi that ¢ appears in combination with the factor 
2U/c, so that we may replace ¢ by the dimensionless 
variable s = (2U’/c)t, which is the distance traveled in 


half-chords. Hence, 


. 1 U —(1/2A)(s—o U 
z= = é w(a) da — W(S) 
Za € FO Ac 


that is, when w 


(9) 


When the gust is a sharp-edged gust 


is constant at the value w,,,,—-we obtain 


, -s/2d 
$= —(UWmar/dAC)e (10) 


The maximum value of the acceleration here occurs at 
the initial instant that the airplane enters the gust and 


is given by 
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GUST LOADS ON 


2 mar. ~ ( W mes. Ac (11) 


Putting this back into dimensional form and dividing 
by the acceleration due to gravity, we obtain the sharp- 
edged gust formula 
(12) 


An pUSW mar (dC,_/da) /2mg 


where ( Am) is the increment of load factor. 

This is, of course, an idealization, and in actuality the 
gust will have a transition zone in which the velocity 
varies continuously from zero to a maximum value. 
The airplane will have attained some velocity in the 
direction of the gust when the region of maximum 
velocity is reached, and there will consequently be some 
alleviation of the gust load when the gust is not sharp- 
edged, even when steady flow is assumed. 

In order to take account of this effect, we now con- 
sider a gust with a transition zone in which the gradient 
is linear (Fig. 1) and js represented by the expressions 


oT 
mar 


lA 


Sg (13) 


Wmar.s S) => Sg 


In the graded portion of the gust the solution becomes 


y ‘ok y 
; 1 l Wmar. (1/2A) (s—o) l W mar 
z= ae da — 
J 0 Ac os 


9\2 ¢ » 
ma COG s 


= (2U/c)(Wmar. / Sq) (e (14) 


For s > s,, the solution may be obtained by subtracting 
from{Eq. (14) the same solution shifted a distance s, 
along the s axis. This solution may be represented 
graphically as shown in Fig. 2. 

It is seen that the maximum acceleration is attained 


when s = s, and is given by 


lad 1) (15) 


Znae * (2U/C)(Wnaz/ 5) ~ 


Dividing this result by the corresponding result for the 
sharp-edged gust {Eq. (11) ], we obtain the following ex- 


pression for the alleviation factor: 


— 970/22 
Alleviation factor = (16) 
$,/2Xr 
The expression for lift on the wing becomes 
l ... 6G * wa) 
L(s) = =pU*S = . —W'(s — a) 
2 da 0 (Cl 


where V(s) is the Kiissner function and ®(s) is the Wag- 
The first two terms are set up by means 
This is 


ner function. 
of the Duhamel, or “superposition,” integral. 
permissible since the aerodynamic theory used is a 
linear theory. The last term represents the “‘apparent 
mass”’ lift and is associated with acceleration of the air 
in the vicinity of the wing.’ Some simplification may 
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gust using steady flow theory. 


The alleviation factor is thus seen to be a function of the 
single parameter s,/A. It is shown plotted in Fig. 3. 

Extending the analysis now to include the effect of 
unsteady flow, we make use of the familiar Wagner 
function,’ which describes the growth of lift on an air- 
foil following a sudden change in angle of attack, and 
the Kiissner function,* which describes the growth of 
lift during the penetration of a sharp-edged gust. 
Both are given in terms of the variable s, introduced 
earlier. They were developed originally for two- 
dimensional! flow but have been modified by Jones to 
apply to rigid wings of finite span.* 


In the present investigation, the functions for two- 
dimensional flow are used in conjunction with the slope 
of the lift curve for the finite wing. While these func- 
tions differ considerably from those for the finite wing, 
their use has been found to yield results that are closer 
to the experimental results obtained from a limited 
number of model tests in a gust tunnel.!” 


i fs : ") d +” pS? 
do + ; (s 7s U do 4 ee 


be introduced by replacing a factor 27 in this term by 


dC, da. 


error and probably helps to take some account of the 


This step, while not rational, introduces little 


finite span effect on the apparent mass lift. 
Transforming to the variable s in all terms, Eq. (17) 


becomes 
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Te ws dC, ‘4 A is ae & where \ is the dimensionless mass parameter defined 
sin 2 ~~ da 0 previously. 
:F° ! - a 
&(s — a)s"(c) do + -2"(s) (18) he exact functions W(s) and #(s) cannot be ex. 
( c pressed in simple analytical form. Satisfactory ap. 
Introducing this expression into Eq. (1) and applying proximations, however, are given by the following ex- 
the Laplace transformation, this time with respect to PFesstons. 
the variable s rather than ¢, but retaining the symbol p 
for the Laplacian operator, V(s) = 1 — 0.500e—-*'!* — 0.500e~° 2] 
a: 1 dC, @(p) P(s) = 1 — 0.165e—%-5" — 0,335e-%39 (29 
ml “( ) p’2(p) = —-= pU*S PrP) ——_ + 
Cc 2 da l 
2 | The Laplace transforms of these functions are 
p?P(p)z(p) + pp) | (19) 
c ( 

0.500 0.500 P 
where a bar over a function denotes its Laplace trans- ¥(p) = p i p + 0.130 - P+ 1 — 
form. This may be rewritten as 

: c [@(p)/ U|V(p) l 0.165 0.335 
2(p) = —- (20) P(p) = —_ —_—_ — (24 
$ p[A + (1/4) + (1/2)(p) ] Pp pt+0.00455 p+ 0.300 


Introducing these expressions into Eq. (20), we finally obtain 2(/P) in the form 


; 0.1412c p* + 0.5756p? + 0.09315p + 0.005141 W(p) 
2p) = - : (25 


A + 0.250 p(p + 0.130)(p + 1)(p? + ap? + aop + a3) Ll 


where 
0.3455 + 0.3364 0.01365 + 0.1438 0.006825 
a, = . a= az = 
+0250 °° d + 0.250 "+ 0.250 
Specializing to the case of a sharp-edged gust, with @(p) = (1/P)Wmaz., 
; 0.1412 ) p® + 0.5756p2 + 0.09315p + 0.003141 a 
a(p) = — : - (26) 
f A + 0.250 l p?(p + 0.130)(p + 1)(p? + ap? + ap + as) 
Noting that 7 = U*(2/c)*s"(s) and that the initial displacement and velocity are zero, we obtain the following 
expression for the vertical acceleration: 
2 = U(2/c)L—' | ps(p)} 
0.5648 Uwnyma: e J p*? + 0.5756p? + 0.09315p + 0.003141 \ ia 
_ > > (24) 


“+ 0.250 c \(p + 0.130)(p + 1)(p* + ap? + ap + az)S 


where £~' | } denotes the inverse transform. The inverse transformation is carried out by application of Heavi- 
side’s partial fractions expansion to yield a sum of exponential terms with real or complex exponents. Thus, 
0.5648 Uwyar 


i. : (Ae 0.1308 4 An |. Bye™ |. Bue? s 4 Bze? ry (28) 
A + 0.250 Cc 


where ¥;, Y2, and y; are the roots of the third-degree polynomial in the denominator of Eq. (27) and 
\ 0.001653 \ 0.5913 
' ~ 0.002197 — 0.01690a; + 0.130a, — a3) ©) 1—a+a—as; 
34+ (.5756y,.2 + 0.09315 + 0.003141 
as Ys Vi ) — ) (k= 1,2,3) 
(ve + 0.180) (ye + 1) (87K + 2diy, + Ge) 
When two of the roots are complex—that is, y2 = a + i8, y; = a — 18—we have 
0.5648 Uwymar , — ae hind cs ; 
=> — = [Aye 0. 130s + Ave + Bie’ + e* (C, sin Ps + Cy cos Bs) | (29) 
A + 0.250 r 


where C, and C, are the real and imaginary parts, respectively, of the quantity 
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(29) 
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(a + 78)* + 0.5756(a@ + 718)? + 0.09315(a@ + 78) + 0.003141 
(a + 18 + 0.130)(a@ + 78 + 1)(a + 78 — y;) 


Dividing Eqs. (28) and (29) by the result obtained in Eq. (11) for a sharp-edged gust under the assumption of 


steady flow, we obtain the acceleration in the form of a dimensionless acceleration ratio. 
0.0648A 


(A,e—9-130s 4. Ave * + Bie? + Boe™ Be" 
+ 0.250 °' ao i 


Acceleration ratio = 
or (30) 


0.5648 " ‘ = er . 
=> [Aje 0.1308 + Ave + Bye’ + e(C, sin Bs + Cs cos Bs) |] 


A + 0.250 


The maximum value of this function is the alleviation factor. It is seen that with unsteady flow there is an 


alleviation effect even for a sharp-edged gust and that the alleviation factor is again a function of the dimensionless 
parameter i. 

When the gust is given an arbitrary shape, the acceleration response may be determined from the response to a 
sharp-edged gust by the application of Duhamel’s integral. This determination is considerably simplified if the 
gust is given a linear gradient, as shown in Fig. |. The acceleration ratio for the period when the airplane is in the 
graded part of the gust is then given by the expression 


0.56480 A, ' . - 
Acceleration ratio = ~ : (1 — e~% 15) + 4/1 —e *) — (l—e™) — 
(A + 0.250).S, L0.130 "1 
B, me B; ak 
(1 — e”’) - (1 — e”) 
Y2 ¥3 
or (31) 
0.5648X f A, ; B, a 
oa > (1 — e~%!85) 4+ 4.11 —e *°) — (1 — e”*) + 
(A + 0.250)s, 40.130 v1 
. i . : C8 + C 
5 —~ (Cia + C28) sin Bs — (CiB — Cre) cos Bs} + aad ay} 
(a? + B*) (a? + B?) ff 


The subsequent response in the constant portion of the gust may be found in the same way as in the previous case 


for steady flow. 


seen that the unsteady nature of the flow produces con- 
siderable alleviation for gusts with small gradient dis- 
tance but becomes unimportant when the gradient dis- 


A program of computations for this case has been 
carried out for a number of values of the mass parameter 
\ and the dimensionless gradient distance s,. These 
two parameters cannot be combined into a single param- tance becomes large. 
eter, as they were in the steady-flow analysis, but must 





be varied separately. 

Acceleration response curves are shown in Figs. 4 to 9. 
It is seen that, for gusts with fairly long gradient dis- 
tance and for small values of \, the maximum acceler- 
ation is reached at the end of the graded portion of the 
gust. This is of some interest, since it is usual, in de- 








termining gust structure from the measured response of 
an airplane, to assume that the gradient distance is 
equal to the distance traveled until maximum acceler- 
ation is reached. For relatively sharp-edged gusts and 
for large values of A, the assumption would not be 
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justified. 

Curves of the alleviation factor are shown in Figs. 10 
and 11. Fig. 11 indicates that the alleviation factor is 
insensitive to gust gradient distance, particularly at a 

















large values of X. This is a consequence of the lag in S 
the growth of lift. The broken-line curves are for the, — ; 

é Fic. 4 Nondimensional acceleration responses using unsteady 
case when unsteady-flow effects are neglected. It is aerodynamic theory. 
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Gust Loaps At HIGH SUBSONIC SPEEDS 


The response to gusts of an airplane flying at high 
speed may be determined in the same way as in the 
preceding section if appropriate expressions for the 
aerodynamic forces are known. For subsonic flight 
there is probably some range of speed in which the 
simple application of the Prandtl-Glauert correction to 
the slope of the lift curve is sufficient to account for 


compressibility. The corrected slope of the lift curve 


may then be written as 
(“*) l dC, 
= (, 

da/x V1— MM da 


where J is the Mach Number and dC,,/da is the slope 
This may be 


to 


of the lift curve for incompressible flow. 
introduced into the expression for dimensionless mass 


parameter A, which then becomes 


mawvVvi— on 
/ - (de>) 
. S d¢ L 
p(c/ =). 
da 


The curves of Figs. 10 and 11 could then be considered 
to apply equally to the compressible flow case, with A as 
defined above. 

The validity of the Prandtl-Glauert correction as it is 
applied here to the gust problem has not been estab- 
Such verification must await comparison with 
the results of a more exact theory. It is not likely, in 
any case, that the correction would be accurate beyond 
a Mach Number of about 0.6. 

In the higher speed range, the work of Possio and 


lished. 


Dietze in nonstationary aerodynamic theory for sub- 
sonic compressible flow forms a basis upon which a 
more rigorous approach may be formulated. Since this 
work was done primarily as an application to the flutter 


problem, it is unfortunately limited in nature. Results 
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are restricted to the forces and moments on an oscillat- 
ing airfoil in two-dimensional flow. This information 
can, however, be used to determine the growth of lift 
and moment following a sudden change in angle of 
attack by applying the reciprocal relations outlined by 
Garrick." Nothing has yet been published on the 
forces and moments due to a sharp-edged gust or a 
sinusoidal gust pattern for high subsonic flow. 


Because of the incomplete nature of the results, the 
existing formulas and tables are not summarized here. 
It may be mentioned, however, that Frazer'! and Frazer 
and Skan'’ have reviewed Possio’s work and give im- 
proved numerical tables for the forces and moments on 
an oscillating thin airfoil for a Mach Number of 0.7. 
The numerical accuracy of Possio’s method is question- 
able at Mach Numbers greater than 0.7 and for values 
of reduced frequency beyond 1. This restriction on 
reduced frequency, although not of great importance in 
the flutter problem, is serious in the gust problem, since 
transient effects theoretically involve all values of re 
Frazer 


duced frequency from zero to infinity. and 


Skan® present results for reduced frequencies up to 2.5, 
although the results for the higher values of reduced 
frequency in this range are questionable, since they are 
based on Possio’s methods. 

It is shown in reference 12 that the lift per unit length 
on a thin airfoil undergoing a harmonic variation in 
angle of attack, of frequency w, may be expressed in 
the form 


P a7 the 
2nr — U*ce™a 


(34) 


i Z3(k, AM) ig 1Z;(k, M)] 





where & is the amplitude of the angle of attack vari- 
(2U/c)t 
is a dimensionless independent variable and is the dis- 
tance traveled in half-chords, and Z3(k, 7) + 12,(k, A) 
is a complex function of A/ and k, tabulated in reference 
12. 

On the other hand, the lift following a sudden change 
in angle of attack, Aa, may be expressed in the form 


ation, k = we/2U is the reduced frequency, s 


L= 21(p 2)U% Aa®, ( M, s) (35) 
where ®,(.\/, s) is a lift-growth function analogous to the 
Wagner function in incompressible flow. 

By means of the reciprocal relations of reference 10, 
the lift-growth function can be related to the functions 


Z3(k, AJ) and Z,4(.\/, k) as follows: 
2 {°Z(k, M) . 
sin ks ds 
TJ 0 k 
ai Z,(k, M) 
TJ 0 k 


It is apparent that the functions 2;(k, J) 
Z,(k, /) must be tabulated over a large range of values 


@,(M, s) = (36) 


or 


cos ks ds_ (37) 


@,(M,s) =1+ 


and 
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Fic. 12. Lift growth function, &,(s). Fic. 15. Lift growth function, ¥z(s). 


of & in order to obtain accurate values of ®,(.\/, s) that a gust analysis becomes practical once again 
from these formulas. The range 0 to 2.5, for which This applies only to two-dimensional flow, since ; 
data are at present available, has been investigated in satisfactory theory for three-dimensional unstead, 
the present work and has been found insufficient to flow is not vet known. 
yield acceptably accurate results. Further work on the - , ; 
; ee ; rhe theory used is an acoustical theory in which the 
gust problem for this case must consequently await the esi ] ; : 
; : : ; : ; airfoil is assumed to be thin and the local downwash i: 
production of more extensive numerical aerodynamic : . ; : . 
let : small in comparison with the velocity of sound. | 
data. na , ; ty ; 
addition, the flow is assumed isentropic and the fluid i 
; ’ ? nonviscous. An expression for lift may be derived 
Gust LOADS AT SUPERSONIC SPEEDS Rite io aa ot i os sas ; , 
which is similar in form to that for incompressible flow 
The available aerodynamic information for super-  Lift-growth functions analogous to those of Wagner and 
sonic flow is much more extensive and, in fact, con-  Kiissner have been presented by Miles,'* Chang, 
siderably simpler than that for high subsonic flow, so Biot,'® and Heaslet and Lomax.'® 


The function giving the growthof lift coefficient following a sudden change in angle of attack of | rad. isdesignate 
?,(s). It is, of course, applicable only for relatively small changes in angle of attack. Whereas the analogow 
Wagner function gives the change in lift due to circulatory effects only, the function ©,(s) for supersonic flow in 
cludes all effects. It is given by the equations 


$,(s) =4/M, O< s < 2M/(M + 1) 


ty Ve | ee | M? —- 1 S lila - : 2 
=e — 2) cos M — VV 5 a > i sin! V/ (1 — =) " 


II 


(58 
; rat Ne 2M 2M 
[t= ar(a- ) a toe 
2M? s Yo ow+i M-—1 
= 4(M?—-1)"”, s>2M/(M — 1) 
and is seen to be a function of the dimensionless variable s = (2U/c)t. It is shown plotted for several values ( 


Mach Number in Fig. 12. 
The corresponding function for lift coefficient following penetration of a sharp-edged gust is given by 


W,(s) = 2s/M, 0<s< 2M/(M + 1) 


Jf ye — 1)-'* [av M? — (5) . 
= : Me? = ) Coe y — 
x \ in uke 2) 


= 4(V/? — 1) » &s22M/(M — 1) 


and is shown plotted for several values of Mach Number in Fig. 13. 
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Using these functions and assuming two-dimensional 
flow, the lift on an airplane penetrating a sharp-edged 
gust and experiencing vertical motion without pitching 
is given by the expression 


1 ow ee Ff eee 
9 = pU*S U W,(s) 4- . pl 2S / 2"(0)®,(s = a) da 
Zz cJ0 
(40) 


where w is the gust velocity, positive up; zis the vertical 
displacement, positive down; and all other quantities 
are as defined previously in the incompressible flow 
case. 

Substituting this expression into the equation of mo- 
tion for the airplane [Eq. (1)] and transforming all 
terms to the independent variable s, we obtain the 


equation 


Sm, w eT & : 
—— ¢°(s) = ——Fz(5) — 3"(o)P,(s — a) do 
pc l cJ0 
(41) 
Rearranging, this becomes 
2"(s) = f(s) — vf 2"(a),(s — a) do (42) 
0 
where 
f(s) = (w/U)(c/4)1/w)¥z(s 
1 /2u 
m = m/p(c/2)S 


is the so-called airplane relative density. 

Eq. (42) is a nonhomogeneous linear integral equation 
of the second kind in the unknown quantity 2”(s), with 
variable upper limit and with an unsymmetrical kernel. 
The kernel &,¢s — o) is not only unsymmetric but has 
the property ®;(s — ¢) = Owhenao > s. 

Solution of Eq. (42) is accomplished here by the 
following approximate method: It is assumed that the 
unknown solution 2”(s) is described by a set of ordinates 
The ordinates 


2", 22", 23”, ..., Zn”, aS Shown in Fig. 14. 
are spaced at equal intervals A along the s axis. The 
integral equation is replaced by a system of linear alge- 
braic equations obtained by letting s and o take on a 
and o; corresponding to the equally 


The evaluation of the 


set of values s, 
spaced stations selected above 
definite integral associated with each selected station s; 
is carried out by the use of the simple trapezoidal rule. 
For example, the area under a portion of the curve 
corresponding to one interval is approximated by 


*h+1 
; l 
/ A(x) dx = » AIA (h) + A(h+1)] (48) 


In this way we obtain the following set of equations: 


3," = fi oo vyAl(1 2)@,(0)s,” + 
k—1 
x é.(s — s)s;"], (k = 1, 2,3, ...) 


j=1 


(44) 


By giving k successive values 1, 2, 3, etc., in Eq. (44), 
each equation can be solved using the solutions ob- 
tained in the preceding equations. 
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Eq. (44) may be made nondimensional by using an 
acceleration ratio as defined earlier instead of 2” as the 
dependent variable. From Eq. (11) we obtain the 
maximum acceleration due to a sharp-edged gust under 


steady flow conditions as 


—(1/p)(dC,/da)(U/c)w (45) 


“mar 


Transforming to the variable s and _ introducing 
1(\/2 — 1)~ ? as the supersonic value of dC,/da, we 
obtain 

Ww ¢ i 
"= — (46) 


U4Yyu?—1 


this into Eq. (44) and representing the 
the 


Dividing 


acceleration ratio by 6, we have set of equa- 


tions 
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bc = ge — yAl[(1/2)4,(0)6, + 
k 1 


>. P,(s, — s;)6;|, (7 — ee 


j=1 


(47) 


where 


. . 
. . 





Solutions of Eq. (48) have been obtained for Mach 
Numbers of 1.25, 1.50, 1.75, 2.00, 2.50, and 3.00 and for 
a selected set of values of u ranging from 0 to ©. Fig. 
15 illustrates typical results for 1/ = 2.0. 

The maximum value of the acceleration ratio in a 
particular case is the alleviation factor. This factor is 
shown plotted as a function of u for several values of 
Mach Number in Fig. 16 A comparison with Fig. 10 
shows that for corresponding values of uw the alleviation 
due to unsteady flow is considerably less in the super 
sonic case than it is in the low-speed case, and, since the 
value of yu is likely to be much higher for a supersonic 
airplane than it is for a low-speed airplane, there is a 
further decrease in the alleviation. In fact, this analy- 
sis indicates that the alleviation due to unsteady flow 
is not likely, in any case, to be greater than about 5 
When the effect of gust gradient is taken into 
account, it will be still less. 

It should be noted again that the analysis given here 


per cent. 


assumes two-dimensional flow and, as a result, may be 
considerably in error in some cases. This requires 


further investigation. 
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Importance of Extending Nyquist 
Servomechanism Analvsis to Include 
Transient Response 


R. L. JOHNSON®* anv J. B. REAT 
Douglas Aircraft Company, Inc. 


SUMMARY 


This article stresses the importance of extending the Nyquist 
servomechanism analysis to include the transient response and 
points out one method of obtaining the transient response from 
the conventional frequency response data which is usually the 
end point of the Nyquist analysis. 

The basis for the present ‘‘rules of thumb’ (phase and gain 
used when only frequency data are considered is re 
a simple 


margins) 
viewed briefly by comparing the characteristics of 
feedback servomechanism with those of a linear second-order 
system. Itis shown that, although these rules are pointed toward 
obtaining a good transient response by considering the system 
frequency characteristics, they often inadequate. Once 
this is realized, the value of actually obtaining the transient 


are 


response is apparent. 

The method outlined for obtaining the transient response to a 
given input consists of representing the input as a Fourier series, 
operating on the series with the system transfer function, and 
summing the resulting series. Response to a step input is given 
is an example, and the necessary procedure for carrying out the 


inalysis is given and discussed. 


INTRODUCTION 


; iow APPLICATION OF Nyquist frequency response 
techniques to servomechanism design is now well 
established. There are several excellent references! ~* 
covering these methods from their mathematical be- 
ginnings to examples of their application. However, 
in using these tools one should not lose track of the 
fact that the analysis is being carried out in a frequency 
domain, while the device itself must operate satis- 
factorily in the time domain. In other words, while the 
analysis and synthesis of the servomechanism may be 
partially carried out in terms of its steady-state fre- 
quency behavior, final assessment of the operating 
device usually depends upon its transient and steady- 
state behavior with respect to time. For example, in 
the response to a step-function change in the input, 
one is interested in the degree of stability, the time 
required to achieve the output within acceptable limits, 
the accuracy of the steady-state performance, and 
whether or not components of the system have been 
These may be easily determined if the 
behavior of the system is known. 


overloaded. 
transient—i.e., time 
The purposes of this article are to stress the importance 
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of extending the Nyquist servomechanism analysis to 
include the transient response and to point out one 
method of obtaining the transient response from the 
conventional frequency response data used in the 
Nyquist analysis. A brief review of frequency re- 
sponse techniques is considered necessary here in order 
to show more clearly the importance of this extension. 
REVIEW OF PRESENT FREQUENCY RESPONSE 
TECHNIQUES 


Prior to the availability of frequency response tech- 
niques, the accepted method for servomechanism 
analysis was the solution of the differential equations 
governing the system." This automatically gave the 
transient behavior of the system, in most cases, to a 
step input. The importance of investigating the tran- 
sient behavior was not overlooked upon the advent of 
frequency response techniques. However, because 
it was often extremely difficult and laborious to solve 
the system equations and because it was not directly 
a part of the frequency response approach, certain 
rules of thumb were developed as substitutes for use 
with the frequency approach. In general, these rules 
are based on experience and on the following reasoning, 
as shown by considering the block diagram of a simple 
servomechanism with negative feedback (see Fig. 1). 


On Fig. 1, 




















uw = forward leg transfer function 

6; = input signal 

6) = output signal, so that % = u(A; — % 
or 

6/6; = w/ (1 + wp) (1) 
6, 
a | 

Fic. 1. Block diagram of simple negative feedback servo 


mechanism 
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90° 


Fic. 2. Nyquist plot for simple feedback servomechanism. 

Determination of the stability of such a system may 
be made by plotting the function » on the complex 
plane with frequency as a parameter. Nyquist’s 
stability criterion for this simple system requires that 
the point —1,j0 be neither intersected nor encircled. 
In Fig. 2, 
curve (b) would be an unstable one. 

Consideration of Eq. (1) and curve (a) of Fig. 2 
shows that the ratio of the output/input amplitudes, 
Bo}, 4:1, is given by the ratio of the length of two 
vectors, one from the origin to the frequency point on 
the curve and the other from the point (— 1,70) to the 


curve (a) would be a stable system, while 
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same frequency point. It can also be shown that the 
phase angle between output and input is the included 
angle between the two vectors. It is evident that for 
example (a) of Fig. 2 a plot of the amplitude ratio ys 
frequency will have the form shown in Fig. 3, where the 
peak amplitude occurs at we = 1, the frequency of 
closest approach to the point (— 1,70). 

Comparison of these frequency characteristics, Fig, 3, 
for a simple feedback servomechanism with those of 
a linear second-order system, Fig. 4, shows 
The amplitude curves both start 


many 
similar features. 
from 1.0 and reach a resonant peak, after which they 
fall rapidly 
start at zero and have lagging angles that chang 


toward zero. The phase curves both 
rapidly in the region of the resonant peak, after which 
Such 


However, 


they flatten out toward some steady value. 
close similarity does not always exist. 
when it does, it leads to the idea of using the transient 
behavior of a second-order system as indicative of the 
behavior of the actual system. This is easily done 
since the relation between the frequency characteristics 
and the transient response of a linear second-order 
system to a step function is well known. Fig. 5 shows 
the transient response of such a system to a step 
function for the same damping ratios as used for the 
frequency characteristics in Fig. 4. 

In using this analogy one usually picks the desired 
transient behavior and, reference to the 
sponding curves of frequency characteristics, notes the 
maximum amplitude rise, the phase lag, and the fre. 
These last quantities 


by corre- 


quency at which they occur. 
are then used directly with the Nyquist plot as criteria 
for the closeness of approach to the point (—1,j0) and 
the frequency at which it should occur. This method, 
however, is only valid if the response of the actual 
system is close to that of a second-order system. 

In an effort to simplify the above procedure, although 
with the same background in mind, a means was sought 
for obtaining the amplitude rise with frequency without 
the necessity of the computation of the quotient y+ 
(1 + yu). This led to the concept of the so-called 
““M circles.”’* Briefly, an .\/ cirele is a curve on the 
Nyquist plot representing the intersection of the vectors 
wand | + wand yielding a constant value of. the ratio 
of the amplitude of these vectors. Fig. 6 illustrates 
several such circles. 

By simply plotting the vector wu, one immediately 
obtains the maximum amplitude ratio and correspond- 
ing frequency by noting the value of the highest .1/ 
circle to which the » curve is tangent. Fig. 6 illustrates 
such a curve, showing a maximum amplitude rise of 
2.0 occurring at a frequency wr & 1.0. 

A further extension of this idea is to plot rectilinearly 
the log magnitude of the vector yu vs. its phase angle, 
again with frequency asa parameter. Fig. 7 shows such 
a plot, consistent with the curve in Fig. 6. 

Some of the advantages of this type of plot are: 
(a) The effect of a system gain change is easily handled 
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py shifting the ordinate scale so that the gain change 
required to bring the system curve just tangent to the 
selected . curve* is easily determined. (b) The .1/ 
curves of constant amplitude ratio u/(1 + yw) need be 
constructed only once, since their geometry is not a 
function of system gain. (c) The reshaping of the 
frequency characteristic required to yield a more stable 
system or the same stability at higher gain is easily 
estimated. 

Rules of thumb have been developed largely on the 
basis of the above reasoning and experience with servo 


and amplifier design using Nyquist techniques. The 


of constant 




















*On a log magnitude-angle plot, the curves 
amplitude ratio are not circles 
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most important of these are the concepts of phase and 
gain margin. These are, respectively, the system phase 
angle at unity gain and the difference between actual 
gain and unity gain at a phase angle of — 180 Fig. 8 
illustrates these quantities on a polar plot 

Common values are 40° to 60° for the phase margin 
and of 10 to 20 db.t for the gain margin. 


observance of these margins not only serves to keep 


Basically, 


down the magnitude of the resonance peak but allows 


+ A decibel is here defined as 20 logio (amplitude ratio 
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, ‘cushion’ for unknown effects and other inaccuracies 


in the data. 


NEED FOR TRANSIENT RESPONSE 


It has been carefully pointed out that the techniques 


apply specifically to those 


namely, those with unity 


mentioned above only 
systems described by Fig. 1 
feedback or nonfrequency-dependent elements in the 


feedback* which obey the relation 
A,/0 u/(1 +p 
in the former case, and 


4/8; (i/K) [uA /(1 + pK) 


in the latter case 

Consider the case in which there are frequency 
dependent elements in both the forward and feedback 
paths. Fig. 9 illustrates this case. 

Subtraction is usually obtained at point A 
i.e., the phase 


NOTE 
by changing the sign of 8 and adding 
ingle of 8 is shifted by 180 


In this case the equation relating output and input is 


6,/0, = iM (1 + uo) 


It may be shown that the Nyquist stability criterion 
still holds—i.e., that a plot of uf should neither inter- 
sect nor encircle the point —1,j0. However, it is 
evident that the relationship between the vectors from 
the origin and from —1,70 to a point on the curve no 
longer furnishes data on the frequency behavior of the 
output /input ratio because the vector ratio yields ug 
instead of w in the numerator. This means that the 
use of 7 circles, ete., is no longer valid, since the 
maximum value of the output/input ratio may occur 
at a much different frequency than that frequency 
where the curve approaches the point —1,j0 most 
closely. As an example, consider the plot of Fig. 10, 
which is a generalized Nyquist plot for the steering 
stabilization of a small missile. Fig. 10a shows the 
general shape of the curve over the whole frequency 
spectrum while Fig. 10b shows the behavior in the 


region of —1,j0. The system is governed by an 
equation of the form 

A uo 

6 1 + wiPi + we82 + 383 


ind the plot shown is the vector sum of the last three 


lenominator terms. 


“It can be shown that the requirement is that there be fre 
juency-dependent elements in only one of the paths—i.e., they 
nay be present in the feedback if there are none in the forward 
Op 

’ For a more complicated system with multiple feedback loops, 
tean be shown that the Nyquist stability criterion is that the 

t number of encirclements, including the contribution of each 


successive open-loop transfer function, must be zero.® 
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While the Nyquist curve passes closest to 1,70 
at frequency fo, a plot of the ratio of output/input 
amplitudes, Fig. 11, shows that the main resonance 
peak actually occurs at f; 

Further, the radical departure of the shapes of the 
phase and amplitude curves from those representing a 
linear second-order system raises a problem of deciding 
whether the behavior illustrated is acceptable. Clearly, 
what is needed is a knowledge of the transient behavior 
of the system. Such knowledge will give the answer 
directly. 

EXTENSION OF FREQUENCY RESPONSE DATA TO 
INCLUDE TRANSIENT RESPONSI 


Fortunately, manual and machine techniques for 
obtaining the transient behavior of a system to an 
arbitrary input, when its frequency behavior is known 
i.e., the reverse process—have recently 
‘4 In essence, the technique re 


and vice versa 
been developed."! 
ferred to in reference 11 is based on matrix methods 
adapted specifically for use with IBM equipment but 
suitable for manual computation also. The first step 
in obtaining the transient output consists of represent 
ing the system input as a Fourier series. Then, multi 
plication of each frequency component of the input by 
the corresponding transfer function of the system yields 
an output component at that particular frequency 
Summation of these output components then yields 
the output of the system in response to the given 
input—4.e., the transient behavior of the system to that 
input. 

As previously stated, the input of greatest interest is 
usually the step input. The equations governing the 
response to a step input have been developed by con- 
sidering the step to be a long-period square wave. If 
the half period of the square wave is selected so that 
the transient dies out before the next step, the result 
is the behavior of a system to a series of step inputs as 
shown in Fig. 12. 

The resulting equation for the output response to a 


square-wave input is: 


AR,, =o 2 ARms 
i, = eee sin |mw,t + Ome 
2 rr ™ 
where 

w, = fundamental frequency of the square wave, 
1/T; 

AR = amplitude ratio, 6 4, a function of fre 
quency 

¢ = phase angle between 6 and @,, a function of 


frequency 

m = multiple of the fundamental frequency (only 
the odd multiples need be considered in this 
case 

Che fundamental period of the square wave must be 


selected large enough so that the transient has died out 
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in the time 7,,/2. 
preferably should not be too large, or excessive labor is 


However, the fundamental period 


required in computing enough terms to make the series 
converge. One way of selecting the fundamental period 
is to estimate the response time based on knowledge 
of the physical situation. Then multiply the expected 
response time by a safety factor of 1.5 to get a value of 
half the fundamental period. If it happens that this 
results in a fundamental period that is too large, and the 
series convergence is thus too slow, then it often pays 
to continue the convergence rather than repeat the 
calculations with a lower fundamental period. It 
should be pointed out, however, that the fundamental 
frequency should always be chosen small enough to 
include the effect of any important narrow frequency 
band peaks in the transfer function. 


It should also be pointed out that the usual Fourier 
precautions relating to boundary conditions and dis- 
continuities should be observed in applying this tech- 
nique to obtain the transient response from the transfer 
The 
precautions should be considered are when the transfer 
function goes to infinity at zero frequency or to a finite 
value at infinite frequency. 


function. most common cases where’ these 


Application of this technique to the data of Fig. 11 
yields the transient response shown in Fig. 13. 


Study of this time behavior shows that the system 
is well behaved, with a well-damped main mode at a 
frequency near f;. The peak in the transfer function 
near fs has put in the higher frequency disturbance 
noted at ¢; and again at ¢;. The steady-state difference 
in output and input is evident from the transient re- 
sponse, as well as from the value of the transfer function 
at zero frequency. The degree of stability, the speed of 
response, and the response time are easily obtained from 
the transient behavior. For example, in Fig. 13, the 
degree of stability is indicated by the shape of output; 
the speed of response is a function of the time, fs, to 
first achieve the steady-state value; and the response 
time may be defined as the time, t, to damp to, and 
stay within, 5 per cent of the final steady-state value. 
Whether or not the system will be overloaded for a 
given input can also be obtained from the transient 
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response. For example, by investigating the transient 
behavior of the system in response to typical inputs,* 
it is possible to determine whether or not overload 
limits are likely to be exceeded. Component behavior, 
as well as overall system behavior, may be investigated 
for overload in the same manner. As a rule, once the 
overall system transfer function and transient response 

obtained the component transfer 
the corresponding transient behavior of 


have been from 
functions, 
any particular component may be checked to see that 
If they 


are not exceeded, the use of the linear theory in the 


the overload limits have not been exceeded. 
frequency response analysis is valid. If they are 
exceeded, one would, of course, have certain reservations 
concerning the results of the analysis. It might also 
be added that the transient behavior of any one com. 
ponent may be of interest from other than overload 
considerations, just as in the case of the overall system 
behavior. 

The work required to obtain the transient behavior 
from the system transfer function is nominal; in most 
cases an average of 3 or 4 hours are required when using 
a desk calculator and 20 to 30 min. are required when 
using IBM equipment or other Fourier synthesizing 
machines. If the amplitude ratio drops rapidly with 
frequency so that not many terms are required for 
series convergence, the time will be shorter than if the 
amplitude ratio curve has an appreciable magnitude 
at the higher frequencies. 

IBM machines are excellent for this type of work. 
Values of the amplitude ratio and phase angle for multi- 
ples of the chosen fundamental frequency are fed 
directly to a set of cards prepared especially for this 
Other 
chanical and electrical machinest for synthesizing a 


purpose, and the machines do the rest. me- 
Fourier series are equally useful for obtaining the 
transient response from the frequency response char- 
acteristics. 

With a little experience in the types of transients 
yielded by typical overall transfer functions, those that 
would yield undesirable transients may be recognized 
and, unless of specific interest, the transient need not 
actually be obtained. This is of importance in a 
synthesis problem where a number of different systems 
require preliminary or rough evaluation. However, 
it is considered axiomatic that final performance pre- 
dictions and system polishing should be based on the 
actual the 
plots and the frequency behavior of the overall system 


transient behavior. However, Nyquist 
should be obtained as logical steps before obtaining the 
transient response for the following reasons: 

(1) 


where trouble is most likely to occur due to uncertainties 


The Nyquist plots point out the frequencies 


* By application of the techniques mentioned in reference 11 
which include the step input as a special case. 

} For example, a machine designed specifically for this purpose 
has recently been developed at the Instrumentation Laboratory 
at the Massachusetts Institute of Technology. 
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ie., in the region of —1,70. The application of the 
criteria of gain and phase margin in this region still 
has merit, especially if these criteria are based on experi- 
ence in relating the transient response to the transfer 
function. If desired, quantitative investigations of 
the effect of unknowns may be made by assuming 
appropriate ranges for the unknowns and investigating 
the effect on the overall system transfer function and, if 
necessary, on the transient. 

(2) Frequencies at which noise and other dis- 
turbances may give trouble are revealed. 

(3) The Nyquist plots provide a convenient graphi- 
cal means of obtaining the vector 1 + yf from the plots 
of u8 and also of obtaining the uf vector for a multiple- 


loop system. 


CONCLUSIONS 


Experience with the evaluation and design of dy- 
namic systems has shown the importance of extending 
the Nyquist frequency response technique to include 
the transient response. It has also shown that the 
transient response may be easily obtained directly 
from the frequency response by the Fourier method. 
Incomplicated multiple-loopsystems with frequency-de- 
pendent elements in both forward and feedback loops, 
the observance of the transient (and associated steady- 
state) behavior is considered to be the primary criterion 
of performance. This point of view is borne out by the 
increasing use of analog and digital computers to yield 
the transient behavior directly. However, for most 
problems, where such large computers are not available, 
the application of the method described herein offers 
definite advantages over the straight frequency re 
sponse approach, even when the use of a desk calculator 
or slide rule is required. When standard IBM equip- 
ment (or a Fourier synthesizer) is available, the manual 
labor of applying this method is small, and its applica- 
tion to many additional aspects of a given dynamic 
problem is made possible without undue expenditure 
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of time and man-hours. It should also be pointed out 
that the manual application of the technique for ob- 
taining the step-function response directly from the 
transfer function is rather routine mathematically and 
may be accomplished by someone not directly connected 
with the servomechanism problem. 
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SUMMARY 


An analysis is made of the effect of skin buckling on longitu 
dinal stability and the influence of such effect on control forces. 
Buckling alters aeroelastic characteristics by the reduction of 
torsional rigidity in the case of an unswept wing and by the re- 
duction of both bending and torsional rigidity in the case of a 
wing. It can also precipitate flow separation, which 
affects pitching moments. 


swept 
Since buckling commonly occurs 
under the higher load factor conditions, large changes in the 
slope of the lines of stick force versus load factor may be the 
consequence. The flow separation appears to have the pre 
dominant influence on such slope changes. Test data of a current 
fighter airplane are used for illustration. 

Mention is made of the effect of a variable modulus of elasticity 
on the rate of structural deformation and of the effect of bending 


deflection on the camber of a swept-wing section. 


INTRODUCTION 


aa LONGITUDINAL STABILITY and control 
determinations are generally based on aerody- 
namic characteristics of smooth-skin surfaces and on 
invariant structural flexibility 
bases are normally adequate, but under high-speed 


parameters. These 
high-load-factor conditions they may lead to erroneous 
results. The factor considered in this paper is the 
buckling of external skins as it affects flow separation 
the 
The 
analysis made of buckling is associated with test data 
on the Republic F-84, a current fighter airplane.! 
A qualitative approach to the subject is made, in view 
of the limited degree of available information. 

F-S4 flight tests, which form the background for this 
paper, were made on an early model of the airplane for 


(and, consequently, aerodynamic coefficients), 


assumption of deformation being linear with load. 


the purpose of obtaining stick-force measurements in 
steady turns at various values of indicated air speed, 
with and without wing-tip fuel tanks installed. The 
empty tank is an unstable body and thus causes a for- 
ward movement of airplane neutral point. For the 
rigid airplane, the pair of tip tanks results in a 2 per cent 
forward movement. An additional movement, vary- 
ing with speed, is shown by conventional aeroelastic 
procedures. Normally, the tanks would be considered 
as simply effecting a reduced slope of the stick-force 
versus load-factor line, corresponding to the neutral 
point shift. However, the test data showed that this 
simple reduction in slope does not occur in the higher 
speed régime of the airplane. 
Received March 28, 1950. 


*Design Engineer—-Technical Liaison, 


ability 


EPSTEIN* 


At speeds up to 500 m.p.h., the stick-force versy 
load-factor lines for both the airplane without tip tanks 
(hereafter referred to as “‘clean airplane’’) and th 
airplane with empty tip tanks are relatively straight y 
through 6g and reduce gradually in slope from there t 
7g, the slope difference between the two airplane con 
figurations being predictable by common aeroelasti 
methods. At higher speeds, to 550 m.p.h. indicate 
air speed, the stick force-trend for the clean airplan 
remains the same. However, when empty tip tank 
are added, there is a reversal in the direction of th 
stick-force line at 5g and 525 m.p.h., and this effect js 
accentuated at 550 m.p.h. Fig. 1 shows plots of thes 
trends and others described below, all taken for th 
same value of center of gravity. 

Flight tests were made with a forward position 0 
center of gravity and with the empty tip tanks installed 
The initial slope of the stick-force-g line was increased 
in line with normal expectations, but the stick-foree 
slope reversal occurred at the same load factor. The 
difference in initial stick-force slopes for the clean air 
plane and for the airplane with empty tanks could 
readily be determined by computation, but no explana- 
tion was available for the sharp demarcation in slopes 
at the higher load factors. 

Flight tests were also conducted with the tip tanks 


full. 


The fuel center of gravity is forward of the wing 
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Comparison of stick force vs. g for various wing-tip-tank 
and speed conditions for the F-84C airplane. 


Fic. | 


50 





elasti 
than 
tanks 
effect 
lead 
450 0 
stick- 
facto’ 
rever 
start 
clean 
the s 
ing ¢ 
slope 
direc 
Fir 
tank 
stabi 
of th 
fin s 
line 
Inste 
verse 
with 
since 
Tl 
wing 
to th 
emp’ 
sion 
tank 
effec 
than 
seein 
acte! 
tank 


Al 
ther 
acte: 
that 
case 
shea 
both 
com 
stres 
of tl 
lowe 
redu 
tors: 
tank 
ativ: 
load 

It 
in tl 
mag 





100 


Versys 
Pp tanks 
ind the 
ight Uy 
here t 
ne con 
oelasti 
dicate 
irplan 
» tanks 
of th 
ffect is 
f thes 
or the 


ion of 
talled 
reased 
<-force 
The 

nN air: 
could 
dlana- 
slopes 


tanks 
wing 


NK 








ank 





EFFECTS OF STRUCTURAL 
elastic axis, and the stabilizing moment of the fuel more 
than offsets the unstable aerodynamic moment of the 
tanks. It was felt that the elimination of the aeroelastic 
effect of the unstable moment of the empty tank might 
lead to norinal stick-force characteristics. However, at 
450 m.p.h., an appreciable reduction in the slope of the 
stick-force line was experienced with increasing load 
factor, and at 525 and 550 m.p.h. the stick-force line 
reversed its direction at 4g (see Fig. 1.), although it 
started out with a steeper initial slope than for the 
clean airplane. The change in slope was approximately 
the same as for the airplane with empty tanks (indicat- 
ing a similar stability change), but, since the initial 
slope was much higher, the pronounced reversal in 
direction did not appear within the limit of the tests. 

Finally, a horizontal fin was placed at the rear of the 
tank with the end in view of contributing some added 
stability to the empty-tank condition. On the basis 
of the usual stability determinations, the addition of the 
fin should merely have rotated the entire stick-force 
line about the origin, increasing the initial slope. 
Instead, the fin eliminated the stick-force slope re- 
versal, and the stick-force lines for the airplane with and 
without empty tip tanks became similar. The fin has 
since become a standard part of the tank assembly. 

The stable torsional moment of the fin about the 
wing elastic axis is approximately equal numerically 
to the unstable moment of the empty tank, so that the 
empty tank and fin combination gives negligible tor- 
sional moment. The stable net moment of the full 
tank without fin, including mass and aerodynamic 
effects, is approximately 40 per cent greater numerically 
than the unstable moment of the empty tank. Yet the 
seemingly more stable full-tank condition was char- 
acterized by a large stick-force slope change, while the 
tank and fin combination was not. 


ANALYSIS OF FLIGHT DATA 


An analysis of these diverse results revealed that 
there is one point in common which appears to char- 
acterize the onset of the stick-force slope reversal 
that is, buckling of the wing skin. Buckling, in the 
case of a wing, is a result of combined axial stress and 
shear stress due to torsion. The empty tank produces 
both increased axial and shear stress in the skin as 
compared to the clean airplane, the increase in shear 
stress being particularly significant in the outer portion 
Thus, buckling effects should occur at a 
lower load factor. The addition of the fin sharply 
reduces the magnitude of the spanwise distribution of 
torsional moment obtained with the original empty 
tank installed. The shear stresses are brought to rel- 
atively low values, and buckling is delayed to a higher 


of the span. 


load factor. 

In the case of the full tank, negative bending moments 
in the outer portion of the span are of appreciably higher 
magnitude than the positive moments for the empty- 


DEFORMATION 


ON STABILITY dl 
tip-tank condition, and torsional moments are higher 
numerically in the outermost 15 per cent of the span. 
The combination of higher axial and shear stress pro- 
duces a slope change in the stick-force-g line 1g earlier 
than for the airplane with empty tanks. The full-tank 
stability change appears due solely to stress effects in 
the outer portion of the span, since stresses are low in- 
board. 

Specific values of the moment distributions for the 
various flight conditions are not detailed here, since the 
varying stress conditions along the span produce dif- 
ferent degrees of buckling, and a definite change in 
aerodynamic characteristics cannot be associated at 
this time with a specific buckle pattern. The buckle 
patterns for the full- and empty-tank (less fin) condi- 
tions are dissimilar except in the tip region; yet they 
produce the same end result. 

So far, the evidence that buckling causes the stick- 
force reversal is circumstantial. A positive check was 
found through the F-S4E, the current model, which has 
heavier gage wing and horizontal tail skins. The heav- 
ier gage sheet produces both a lower stress level and a 
50 per tent increase in computed buckling stress. 
Thus, under 6g flight conditions, the surfaces remain 
Flights were made at 550 m.p.h. with empty 
In contrast to the 


smooth. 
tip tanks without fins installed. 
stick-force slope reversal at 5g which had been pro- 
duced by the thinner skinned airplanes, the stick-force 
line for the F-84E, under identical flight conditions, was 
found to be relatively straight up to 6g, the highest load 
factor used in the tests. 

Buckling can produce two effects 
rigidity and the onset of flow separation at a lower Mach 
Number than that of a smooth surface. The reduced 
torsional rigidity leads to greater positive wing twist, 
which in turn increases wing slope and reduces longi- 
Flow separation at the critical Mach 


reduced torsional 


tudinal stability. 
Number of the F-84 wing induces a nose-up pitching 
moment. Earlier flow separation caused by buckling 
could induce this nose-up effect at a corresponding 
Mach Number. The stick-force slope reversal could 
then follow, since a movement of the elevator in the 
down direction is required to balance the nose-up mo- 
ment. One clue in this regard is the behavior of an 
F-84 having a few leading-edge panels slightly dished 
in. The pitching-up tendencies characteristic at critical 
Mach Number occurred 0.03. earlier. 

Adams and Silsby? have shown that a faired bulge 
with a maximum depth of 0.13 per cent of the chord, 
extending from 35 to 55 per cent of the chord of the 
wing of one of the wartime fighter airplanes, reduced 
the critical Mach Number by 0.01. and had a great 
effect on pressure distribution when the critical Mach 
Number was exceeded. Rhode* compared the incom- 
pressible pressure distributions of a smooth wing chord 
section and of the corresponding wrinkled section 
formed in static test. Appreciable increases in the 
maximum negative pressure coeflicients are shown to 
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result from the wrinkling. In the case of the empty 
tank condition of the original F-S4 model, the aero 
dynamic effect of the skin wrinkling can be ascertained 
by comparing the known critical Mach Number of the 
smooth-skin wing with the Mach Number where a large 


change in stick-force slope occurs at high g. Flow 
separation is indicated at least 0.05.7 earlier. If de 
parture from linearity of the stick-force line at 450 


in.p.h. for the full-tank condition is taken as the criter 
ion, then the buckling effect on aerodynamic character 
istics extends downward to a value 0.15.17 below the 
critical Mach Number for the smooth surface. 
sibly, the more pronounced buckle pattern for the full 


Pos 


tank condition is the cause of the earlier aerodynamic 
transition. 

An estimate of the effect of reduced torsional rigidity 
can be made. The remaining influence on the stick 
force slope reversal may then be attributed to the earlier 
induced flow separation In the case of the full-tank 
tests, the net torsional moment about the wing elastic 
axis is relatively small, except in the outer portion of the 
span where it is stable, and thus rigidity changes cannot 
be a factor for the full-tank stick-force slope reversal. 
Since the extent of the influence of such changes is 
desired, only the empty-tank condition is considered in 
further discussion. 

Two investigations of the N.A.C.A.*° yield data 
showing the effect of buckling on torsional rigidity. 
Reference 4 reports twist versus torsional moment for a 
box of section 24 by 10 in., having 0.026-in. skin, a 
longitudinal stiffener spacing of 4 in., and a transverse 
stiffener spacing of 19 in. The onset of buckling as 
determined by strain gages was at a stress of 2,660 Ibs. 
per sq.in., as compared to 2,260 and 3,790 Ibs. per 
sq.in. for simple and clamped supports, respectively. 
At a torsional moment double the value causing the 
initial buckling stress, the loss in rigidity was 30 per 
cent. At still higher torques, further reduction in 
rigidity was obtained. 

Reference 5 is concerned with cylinders with various 
degrees of skin stiffening. As an example, for a 30-in. 
diameter cylinder having 0.020 24S-T skin with 9-in. 
peripheral stiffener spacing and 5.9-in. longitudinal 
stiffener spacing, the loss in rigidity after buckling 
This is equivalent to 
In the 
case of a wing, the contour behind the usual location of 


was approximately 75 per cent. 
a deflection rate four times the initial value. 


front shear web has a slight curvature; the contour of 
the leading edge had a considerably greater curvature. 
A box uniform along its length, using a typical wing 
section, would experience a greater loss in rigidity after 
buckling than the flat-sided box, but it is not believed 
that the difference would be pronounced. 

A wing is a more complex structure than the speci- 
mens described above. The unit skin-panel dimen 
sions, which control] the magnitude of critical buckling 
stress, vary along the span, and the stress condition, as 


already noted, is one of varying combined axial and 
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shear stress rather than of simple shear from a torsional 
moment. An analytical determination of the varia 
tion of torsional rigidity with load factor thus appears 
unpractical. Because of the minimum gages of skip 
required for reasonable aileron power, flutter, ete., , 
modern fighter wing does not normally buckle in th, 
lower load-factor region, and thus the use of constan; 
deflection rates in aeroelasticity is reasonable in mog 


applications. 


EFFECT OF RIGIDITY CHANGE 


The static test of one of the early models of the F-s 
airplane furnishes an example as to the loss in rigidity 
that might occur in an actual wing as the load or loa 
factor is taken it 
merements of 20 per cent of ultimate load, show n 


increased. Twist measurements, 
significant change in rigidity up through limit load 
The torsional deflections in the increment from 60 t 
80 per cent of ultimate load run approximately 20 pe 
cent more than in the increment from 20 to 40 per cent 
and the deflection rate in the increment from SO to 9 
per cent of ultimate load runs approximately double th 
average rate from 20 to SO per cent. It should be note 
here that the technique usually used in obtaining dat; 
applicable for determining differences in rigidity over 
small load increments is of limited precision, and, ther 
fore, the loss of rigidity should be regarded more as 
trend than as quantitative data. 

This particular static test was conducted by loading 
shot bags on the wing. Small shifts in placing the bag 
from one increment to the next and limitations in read 
ing scales may introduce considerable error in the tor 
sional deflection results On the other hand, loading 
through tension pads (steel-backed rubber pads ce 
mented to the wing surface) may introduce error through 
pad restriction on free-buckle formation. 

The static test condition was for the clean airplane 
When the empty tip tank is added, the bending an 


torsional moments are increased along the entire span j 


the greatest percentage increase being over the outboar 
section. The increase in net bending moment at th 
mid-semispan section is such that the bending moment 
at this section at limit load is of the same order of mag 
nitude as it is at 80 per cent ultimate load for the cleat 
wing. 
the moment of the clean wing at the mid-semispat 
section at SO per cent ultimate load is realized at les 
than limit load factor. Thus, at limit load, this sectiot 
would be under a stress condition that causes a sub 
stantial loss in rigidity, and sections further outboar 
would have the rigidity loss at less than limit load. 

A precise reduction in rigidity for the loading condi 


tion with empty tip tanks installed at, say, 5g cannot 


be determined on the basis of available data. At! 
indirect approach was therefore used in determining t 
what degree nonlinear deformation rates can be hele 


accountable for the stick-force slope reversal at 59! 


The increase in torsional moment is such that] 
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mp h. The following describes the procedure followed 
Since previous discussion has shown that the reduction 
in rigidity 1s likely to be over the outer portion of the 
span, it was assumed that a constant degree of rigidity 
loss applies over the outer half of the semispan and no 
loss over the inboard half. 

First, the spanwise variation of rigidity determined 
froma torsional moment applied at the tip was assumed. 


The 


tip moment is generally too small to produce any sig 


[his is the rigidity normally attributed to a wing. 


nificant buckle pattern except in the immediate vicinity 
of the tip 

The second assumption was that the outer half of the 
panel experiences a 30 per cent loss in rigidity, and the 
final assumption was that the outer half experiences a 
50 per cent loss in rigidity. These are arbitrary figures 
selected to show the significance of a large rigidity loss 
on longitudinal stability. Even the 30 per cenit loss 
is considered high for the 5g case. 

For the first rigidity condition, the wing twist was 
change in root angle of attack at 
The Multhopp span 


determined for a ] 
550 m.p.h. indicated air speed. 
load procedure was used in conjunction with a process 
The resulting total lift on the span divided 
angle-of-attack change 


of iteration. 
by the rigid wing lift for a 1 
gives the increase in lift-curve slope due to the twist. 
[he same procedure was followed using the other two 
rigidity assumptions. Increased lift-curve slopes were 
obtained, and comparison with the slope for the first 
rigidity case gave the percentage increase in slope to use 
in obtaining neutral-point changes due to assumed 
rigidity losses 

The first lift-curve slope is the one responsible (in 
conjunction with corresponding fuselage and tail de 
formations and aerodynamic parameters) for the initial 
slope of the stick-force line for 550 m.p.h. for the empty- 
tip-tank case. For a 30 per cent loss in outer panel 
rigidity, the wing lift slope is increased by 7 per cent, 
and for the 50 per cent loss, the wing lift slope is in 
creased by Using first-order effects on 
stability, the neutral-point shifts amount to roughly | 
and 3 per cent of the M.A.C. for the 30 and 50 per cent 
To account for the full 


l7.0 per cent. 


losses in rigidity, respectively. 
change in stick-force versus g slope from the initial value 
to the average value between 4.5 and 5.5 g would require 
a neutral-point change of 7.5 per cent. Thus, even an 
extreme loss in torsional rigidity can account for less 
than half the overall effect of buckling of the F-S4 con 
trol forces. . 

The balance would then logically be the equivalent 
of a change in wing moment coefficient due to the early 
flow separation A coefficient change of 0.005 is all 
that is required to cause the complete stick-force slope 
to 5.5g load 


the 


reversal effect at 550 m.p.h. in the 4.5 


factor increment \t lower indicated speeds at 
Number, the effect 


since less elevator deflection is required to offset the 


same Mach would be less severe, 


same moment coefficient change, and for a particular 


DE 


FORMAT ON STA 


LON 
elevator deflection the stick forces are reduced in pro 
portion The for lowet 
elevator deflection follows from the smaller reduction o/ 
to effects at the 


to dynamic pressure. need 


effectiveness due aeroelastic lower 
speed. 

Factors that remain to be settled regarding the early 
flow separation are the quantitative effects on moment 
coefficient, lift 
cific buckle pattern, and the relation between depth 
of buckle the the 


cients. 


coefficient and downwash of a spe 


and severity of change of coeffi- 


OTHER STRUCTURAL DEFORMATION EFFECTS 


The nonlinear rigidity factor resulting from buckling is 
not limited to the wing. It could also apply to the 
horizontal tail and the fuselage. 
may be a factor not only in longitudinal stability and 
control but The 
magnitude of the rolling power of a displaced aileron 
may be influenced by skin buckling. The early flow 
separation aspect of buckling can apply to the horizontal 


Also, the nonlinearity 


in other phases of aeroelasticity. 


tail with consequent change in elevator effectiveness and 
hinge moment. 

In the swept-wing airplane, bending of the wing be 
comes a factor in the change in neutral point, since it 
affects the section angle of attack in the direction of the 
air stream. Thus, buckling of the skin at high load 
factors could reduce the bending rigidity of the wing 
and increase the unstable moment introduced by wing 
bending. Of course, premature flow separation can 
affect a swept wing as well as a nonswept wing. 

A factor involving structural deformation, which is 
peculiar to swept wings, is the change in section camber 
caused by bending deflection. This change in camber 
results from the change in slope of the wing along its 
length from the leading edge to the trailing edge of a 
chord line taken in the direction of the air stream. If 
stations are measured along the wing 50 per cent chord 
line rather than along a line normal to the plane of 
symmetry, the leading edge may be at, say, station 100 
and the trailing edge at station 140. The offset 
tween a line joining these two stations and the curved 


be 


line of the deflected wing between these stations gives 
the camber offset, and the distance divided by the chord 
length gives the degree of camber of the new airfoil 
The usual aeroelasti¢e treatment of the swept wing con 
siders the basic airfoil shape of the rigid wing as remain 
ing unchanged, even though it no longer exists after the 
wing is bent. A check of the camber change for one 
swept-wing fighter airplane shows that 
one-half of | per cent is attained at limit load factor 


a change of 


over the major portion of the span. This camber figure 


is not necessarily a typical one. At a particular load 
factor, the camber change will be different at supersonic 
alt 


speed from that at subsomic speed because the 


movement of aerodynamic center results in a different 


twist condition 
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In reference 6 for the 4-digit-series N.A.C.A. airfoils, 
a change in camber of | per cent, positively, is shown to 
change the zero lift angle by plus 1° and the moment 
coefficient by —0.025, at a low Mach Number. Since 
positive bending produces negative camber, the effect 
of the bending is to raise the value of zero lift angle 
(which reduces the wing lift curve slope and is stabiliz- 
ing) and to induce a stalling moment of increasing 
magnitude as the wing bends upward. Both of these 
factors would modify the deflected shape of the wing 
as determined by conventional aeroelastic procedures. 
The effects are gradual and would not produce a sharp 
change in stick-force slope shape as is caused by buck- 
ling. Experimental verification of the influence of 
bending on aerodynamic parameters through camber 
change would be of interest. 

Mention should be made of another nonlinear aspect 
of structural deformation—the variation of the modulus 
of elasticity with stress. In the case of Alclad sheet, 
there is a drop of 10 per cent in the value of the modulus 
at 10,000 Ibs. per sq.in. tensile stress. The proportional 
limit of 24S-T sheet is less than half the ultimate tensile 
strength, and a still greater disparity exists for 1S8-S 
corrosion-resistant steel. However, for the usual stress 
levels within permissible flight conditions and for the 
normal proportions of clad to unclad material in a wing, 
the modulus change is considered only an academic 
point. 


CONCLUSION 
This paper has discussed some aspects of structyr 
deformation which appear to deserve serious consider; 
tion in longitudinal stability and control estimate 


Even if skins are designed to be nonbuckling up to lim; 


load factor, some concern should be given to stick-forg 
changes that might occur if the limit load factor ; 
exceeded inadvertently The treatment of the subje 
has perforce been sketchy. Additional research coul 
fill some of the obvious gaps in the data given in th 
paper. 
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Steady, One-Dimensional, Viscous and 


Compressible Gas’ 
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SUMMARY 


In this paper the equations of motion for a steady, one-dimen- 
sional, viscous and compressible gas are simplified according 
to two procedures and are then integrated. In one procedure the 
viscosity term is retained in the momentum equation and omitted 
in the energy equation. In the second procedure the viscosity 
term is omitted in the momentum equation and retained in the 
energy equation. 

The equations used in the first case are the momentum equa 
tion, into which is introduced the frictional force [= (4/3)ud?u- 
ix*|; the continuity equation; and a polytropic relationship 
between pressure and density, the last expression replacing the 
energy equation in its usual form. The value of the exponent n, 
in the polytropic relationsbip for which the Rankine-Hugoniot 
points are satisfied, is shown to be a function of the initial Mach 
Number, 1/9, only. It is strongly indicated that retaining the 
viscosity terms in the momentum equation and neglecting them 
in the energy equaticn would offer a significant and considerably 
simplified mathematical representation of a steady compressible 
viscous flow in two and three dimensions. 

In the second case, the Eulerian equation, together with the 
continuity and energy equations for a one-dimensional motion, 
describes approximately the flow of a viscous, heat-conducting 
gas. This is analogous to the treatment of Boley and Lieber! 
for two-dimensional flow. 

An exact solution is obtained for the shock-wave structure in 
both cases, and, although the velocity distributicn in the latter 
case is continuous everywhere, the singularity in the velocity 
gradient is not removable. However, while the present report 
shows that retaining the viscous term in the energy equation 
ilone is not sufficient for the removal of this singularity, it is shown 
that the frictional force, when considered in the momentum 
equation, does remove the singularity. Since the shock-wave 
thickness, according to the Prandtl definition, is zero in the second 
treatment (a physical impossibility), it was not found necessary 
tocalculate the variation of the remaining flow variables (p, p, S, 
l) with Mach Number. 

Good agreement with the exact one-dimensional solutions of 
Morduchow and Libby? and H. Reissner and Meyerhoff® with re- 
gard to the structure of the shock wave and its thickness is ob 
tained when the effect of viscosity is considered only in the mo 
mentum equation. Moreover, indications are such as to justify 
the conclusion that the presence of the viscous stress term in the 
momentum equation is a necessary condition for the removal of 


the singular velocity gradient in one-dimensional flow 
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INTRODUCTION 


§ ie PURPOSE OF THIS PAPER is to evaluate two 
approximate methods for representing the effects 
of viscosity in the equations of motion for a compressible 
viscous flow. These investigations are carried out 
here in one-dimension, since the solution to the equa 
tions of motion is readily obtainable in this case and 
can be compared to available solutions of the exact 
equations of motion. Consequently, the solution ob 
tained from the simplified equations of motion con 
sidered here in one-dimension is not intended for the 
purpose of obtaining new information but rather of 
establishing a basis for making corresponding simpli 
fications in two and three dimensions where the exact 
equations of motion are, in general, mathematically 
intractable. 


In a previous report by Boley and Lieber! on two 
dimensional shock waves, the viscosity terms in the 
momentum equation were neglected. Consideration 
of energy dissipation by viscosity was shown to be a 
sufficient condition for the elimination of the singulari 
ties arising in perfect fluid solutions in the neighborhood 
of the local speed of sound. The results of that in- 
vestigation are valid only for Mach Numbers approxt- 
mately equal to unity. This is due to the assumption 
that the perturbations of the velocity components are 
small, though their gradients in the direction of flow 
are large. 

The position and direction of the shock were obtained, 
although the structure of the wave is not determined 
by such an analysis. The reason for this is probably 
the fact that the viscosity terms in the energy equation 
appear as scalar quantities. The treatment in refer- 
ence | requires further elucidation, and the results are 
tentative. 

Since the frictional force introduced into the momen- 
tum equation is of vector form, it was expected and is 
shown as true that detail on the nature and structure 
of the shock wave could thus be obtained. 

It is therefore believed that retaining the viscosity 
terms in the Navier-Stokes equation and neglecting 
them in the energy equation would provide significant 
mathematical descriptions of flows in two and three 


dimensions as well. 
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l subscript denoting values at ¢ 0, inflection 
es pom : Substituting Eq. (6) into Eq. (5) gives 
Cy, Ce, Cs, Ce constants of integration 
C. = specific heat at constant pressure | du du 
C. specific heat at constant volume punt! L (Anm mu" +!) 0 ’ 
a 
R gas constant = C, c 2 dy an whe 
1 = ratio of specific heats C,/C ; : 
h (euitiant ef tant coniuctivits and applying the transformations q¢ du dx and q(dq Ir 
V Mach Number du) d°u dx* to Eq. (7) gives L.e., 
n exponent in pressure-density relationship, 
p = Ap" $ du Am : 
a a mu C; \ 
p pressure os @& u! prov 
p mass density of fluid 
> te » . ‘ cons 
l absolute temperature Eq. (S) can now be integrated into the form 
rv coefficient of viscosity Eq 
u velocity of flow in x direction | > ude . 
| dimensionless velocity u/% \ L. ¢ 
Ss entropy 3 J mu"t!+ Cu" + Am’ tt 
° . as 
thickness of shock wave, according to Prandtl ; 
definition It was found convement to nondimensionalize bot It 
2 distance in direction of flow wu and x with respect to up and Xo, respectively, wher ft 
4 . . . . 0 
r = mean free path of gas molecuk the subscript zero denotes values at minus infinit —_ 
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[2/(y + 1)]M 
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) : { ( \ vk 
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introduced only into the momentum equation, are: ; 
' | No { Un Z po) Vir SR 0 Fi 
Conservation of momentum: 7 ; ‘ , , ber, : 
' Using relations (10) and (11b), the integral (9) can ’ 
) ) . . . : the \ 
du dp $ du put into the nondimensional form 
ui — = - j 7-2 ly? (1) IS ap 
. x 3 x r 
ax ax ) ( | IsR7 c Macl 
. . . . <¢ = Z po p 4 
Continuity equation: 3 s 5 may 
d dx)(pu) = O (2) [ up" t1V" dV 
. i P 7 : n+l n+) r ’ ni , 
Pressure-density relation: SJVi MU + Cyuo"V" + Am 
wher 
p = Ap’ (3) where 
In the above equations, u is the velocity in the x- A Po) po" = RT / po’ 
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Rankine-Hugoniot points are satisfied. 
; (9) ov = ee See ee a . a z 
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pu = pollo = m (4) ing : arv c ‘ > slope 
p Apply ing the bi undary condition th it the sk P given 
ver = : the velocity distribution is zero at minus infin 
Substitution of Eq. (4) into Eq. (1) leads to ; sett, : 
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{ SZ l al TT 
= u 
$ ; yr Mo St“ 2 ioe,’ ++ +———+ 4 
“V 
d\ | | 
/ 4 (14) +\ + 4 
ei : 
(| 1) : (1 z -+——+- | } + 4 | | 4 
yMo y" | 
= Se - Oa a 
where W/o uo/ V yRT >. 
In order to locate the origin at the inflection point 
ie., du/dx? = 0 when x = 0) Eq. (7) yields | — | ‘2 ice a os ee | 
uw, = (Anm"—')/% +0 (15a) t — a on an a a | 
{ | 06 
provided « # constant. This is valid because wu r i oe oo | [ Ls | T 
constant is a trivial solution. Nondimensionalizing L [ —- - } ss | } 
; ay rives | | 
Eq. (15a) gives aa =— te 5 | wiscosity TERM—1 
, ry 1/(n + 1 : RETAINED IN THE 
Vi (n/y Mp?) (15b) | MOMENTUM EQUATION | 
im | | ONLY 
- ° ° . P | 
as the nondimensional velocity 1); at & = 0. a = Leal | | , 4 fs 
: | 
It is obvious from Eq. (14) that IV’ | satisfies one q | | | hee in, bet ete 
| 
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of the Rankine-Hugoniot points at minus infinity 


In order to satisfy the second point, 


, es 2 
| a t 2 
+ I (y + 1) Mo 

the value of ” 1s found for which the denominator of 

Eq. (14) will vanish. Hence, 
y¥M?(a 1) | (1/a") (16a) 

Solution for n gives 

n = —}log [1 + y\W02(1 a)| log af (16b) 


Fig. 4, which is a plot of m versus initial Mach Num- 
ber, shows that » > y for all Mach Numbers. 


the value of the Mach Number at the inflection point 


Hence, 


is approximately equal to ie eo the initial 


Mach Number is not too large M, < 3.0. This 
may be shown as follows: 
M=u/VyRI = VM,/V T/T, (17a) 


uyV and M = wV yRTp. 


T/T» = 


where u However, 


pu/poto = 1/V"—! 


and substitution of this expression into Eq. (17a) leads 


to 


M = MV V"t! (17b) 


The value of the Mach Number at the inflection point 


isa function of the nondimensional velocity, V, which is 


given by expression (15b). Hence, 
M, = Vn/¥ (17c) 
Thus, the Mach Number at the inflection point, 
M,, is approximately equal to unity when n ~ y—i.e., 


for M, < 3.0. 


Figs. 1 and 2 show the velocity distribution for M 
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FIG 2 VELOCITY DISTRIBUTION 


bers with those of Morduchow and Libby® shows good 


agreement. 


Shock-Wave Thickness 


The thickness of the shock wave (that region in which 


the major portion of the change in the flow variables 


occurs) as defined by Prandtl is, mathematically, 


Uy ee 


(du/dx) max 


(1S 


Noli 


Nondimensionalized, Eq. (1S) becomes 
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FIG. 3 SHOCK WAVE THICKNESS BASED ON 
PRANDTL DEFINITION 
th = (a — 1) (dV /dé), - (19) 
where 


* 
dé t=0 
3 lyn al . | ( -) | 
Vy,-1) - lL — 20) 
4 s PA ; 7M? Vi" ( 


Here, VV; is defined by expression (15b). Fig. 3 
shows the comparison for shock-wave thickness, based 
on Prandtl’s definition, with those of Morduchow and 


Libby.” 


Entropy 


It will now be of interest to determine the change in 
entropy, AS, asa function of V. In general,” 


- T p 
= C, log — — R log (21a) 


T» Po 


Using the relations C, = yC, and R = C,(y — 1), to- 
gether with the expressions p/po = (p/po)” and 7/7) = 


(p/po)"—', Eq. (21a) becomes 


Ss >= So = AS 


S — So 


- = (y — n) log (”) = (y—n)log V- (21b) 
Ci p 


where n is determined by Eq. (16b). 
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FIG. 4 VALUES OF n NECESSARY TO SATISFY 
RANKINE - HUGONIOT CONDITIONS 


Eqs. (21b) and (14) relate the change in entropy 
(AS/C,) to — as shown in Fig. 5. A plot of m versus 
initial Mach Number, Mo, in Fig. 4 shows that m > 
for My = 1.0. Hence, the change in entropy is always 
positive in accordance with the second law of thermo- 


dynamics. 


Viscosity EFFECTS INCLUDED IN THE ENERGY 
EQUATION AND OMITTED IN THE MOMENTUM 
EQUATION 


The Eulerian equation in one-dimensional form, to- 
gether with the continuity and energy equations, may 
be used to describe a pseudoviscous heat-conducting 
gas. These are: 

Momentum equation: 


pu(du/dx) = —(dp/dx) (22 
Continuity equation: 

(d/dx)(pu) = 0 (23 

Energy equation: 

ee du 4 du\? dT 

ee dx TP dx 7 3 “ (“*) ~* dx? “eo 

Equation of state: 
b = pRT = mRT/u (25 
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Here the flow variables, p, p, 7, and u are defined as Integrating Eq. (26), one obtains 
before. C, is the specific heat at constant volume, R > + mu = Cy = po + its (28) 
isa gas constant = C, — C,, and k is the coefficient of 
heat conductivity. Eq. (28) can now be used to eliminate p from Eq. 
Integrating Eq. (23) and substituting this into Eqs. (27). From Eq. (25), an expression for T and, hence, 
(22) and (24) gives T’ and 7” can be obtained, where the primes denote 
differentiation with respect tox. Thus 
m(du/dx) = —(dp/dx) (26) 
T = (1/mR) (Cyu — mu?) 
C dT 1 ,du 4 (“ty 7 pel =0 (97) T’ = (1/mR) (Cyu’ — 2muu’) (29) 
wel” f dx 3° \dx dx? - T” = (1/mR) (Cyu” — 2mu’? — 2muu") 


Substituting the expressions for p, 7’, and 7” into Eq. (27) gives, upon nondimensionalizing uw and x, 
C, (: ‘i ) . ve lt() 4 l |e 4 (~) 
R B2/ dé —— dé B2 dé  3(R.N.) \dé 
C, 1\ @V dV \? €V 
i+ ) - 2 ) — 2] ~1=0 (30) 
. RP,(R.N.) B2/ dé dé dé? 


B = Y yMo, R.N. = MX fo = Vv yr 8(Mo, SAN P. = uC,/k 


where 


Letting g = dV /dé and q(dq/dV) = d?V/dé?, Eq. (30) becomes 


(dq/dV) + g[C/(aV + b)| = (eV + f)/(eV + d) (31a) 
where 
a = (2C,/R)/PA(R.N.); 6 = —al[l + (1/82)]/2 
4 1+ Cc, , (: Fi “)( ‘ *) (31b) 
, =a ‘ e= : = — 
3(R.N.) g* B? Rr) 


Eq. (31a) is a first-order linear differential equation that can be integrated to give 
q = [e(aV + 6)/a(a + c)| + (l/c) (f — eb) + [Cy/(aV + 6)°*] (32) 


where C, is determined by the boundary condition 


g=OatV=1 (33a) 
Hence, 
Cy, = —[e(a + b)' t+ “® /a(a + )| — (1/c)(f — eb) (a + 0)" (33b) 
Thus Eq. (32) becomes 
e(aV + bd) l a+b\““*[e(a + b) l 
=” = 2 * c laa. SS ;) Ee : Cc) . c Y — )| — 


It can be seen from Eq. (34), since c and a are both positive and b is negative, that —b/a < V < 1. Since 


“Vy 
&= | dv/q 
Vi 


then VY; = —b/a. Thus, the velocity distribution is defined from the origin to negative infinity. In order to 
obtain the velocity distribution from the origin to plus infinity, the other value of the Rankine-Hugoniot point is 
utilized as a second boundary condition, 


gq=OatV=a (35) 


Hence, Eq. (32) becomes 


e(aV + b) rt l f i (2 4 a Ee + d) 1 | (y ; | 26 
= . iTt-  =— — €0) (36) 
q a(a + c) "aad P aV+b a(a + c) eC 


wherea < V < Vj. 
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FIG 6 VELOCITY DISTRIBUTION 


In order to obtain the velocity distribution as a 


function of &, the expression 


is used, where |’; = —)d/a. 

According to the Prandtl definition of shock wave 
thickness given by Eq. (19), a zero thickness is ob- 
tained since (dV/d& 29) = — ©. This physically 
impossible condition, shown in Fig. 6, is sufficient reason 
to discard any calculations regarding the distribution of 
the flow variables through the shock. However, the 
velocity distribution is continuous, showing that the 
discontinuity present in perfect-fluid flow analysis may 
be removed by the consideration of a frictional term 


only in the energy equation. 
CONCLUSIONS 


(a) Retaining the viscosity term either in the mo- 
mentum equation or the energy equation is sufficient 
to remove the discontinuities in the flow variables 
which exist in the analysis of a potential flow. The 
latter case, however, because of the presence of an in 


finite velocity gradient, is not to be accepted. 


60) 


(b) The velocity distributions, with viscosity ¢op 
sidered only in the momentum equation, are in gog 
agreement with those of Morduchow and Libby? an 
H. Reissner and Meyerhoff,* despite the fact that th 
Prandtl Number assumed in this report is infinite 
i.e., Rk = (-—as compared to a value of */, used in th 
above references. Hence, it is indicated that Prand 
Number has little effect upon the structure of the shog 
wave. 

(c) The agreement in (b) indicates that retainin; 
viscosity terms in the Navier-Stokes equations an 
neglecting them in the energy equation may give 
simplified, although physically significant mathemat; 
cal, description of two- and three-dimensional, con 
pressible, viscous flows. 

(d) The value of mn, exponent in the pressur 
density relationship, is shown to be a function of th 
initial Mach Number only if the velocity distribution 
to satisfy the Rankine-Hugoniot points. Since n >, 
for JJ > 1.0, the second law of thermodynamics is satis 
fied—1.e., compression waves are accompanied by pos 
tive increases in entropy. 

(e) The shock-wave thickness as shown in Fig 
varies inversely with the initial Mach Number andi 
expressed in terms of a pure number—namely, th 
number of mean free paths at minus infinity. Sin 
there is little doubt that the continuum theory wou 
be valid if the thickness were large, it can be stated that 
for Mo < 1.25 and¢ > 6.0, the equations used to describ 
the flow are adequate. 

(f) The infinite velocity gradient cannot be remove 
when the viscous stress term in the momentum equa 
tion is neglected. This was shown to be true wher 
first the viscous and then the heat-conducting tem 
in the energy equation were omitted separately. Neg 
lecting the viscous term and retaining the heat cot 
duction term in the energy equation resulted in a solu 
tion similar to the one obtained in Fig. 6—1.e., a singu 
lar velocity gradient. When the heat-conductio 
term was neglected and the viscosity term retain 
only in the energy equation, the differential equatio 
was reduced to first order. Hence, the solution cou 
not satisfy the boundary conditions. 
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An Approximate Method for the Calculation 
of the Stresses in Sweptback Wings 


THORKILD RAND* 
Royal Institute of Technology, Stockholm, Sweden 


SUMMARY 


An approximate method is proposed to eliminate the error 
that arises when groups of stringers are replaced by a few imagi 
nary longerons in order to reduce the number of statically inde 
terminate quantities in monocoque structures. The idea under 
lying the method is to choose arbitrarily a limited number of 
degrees of freedom of motion and to consider the corresponding 
self-equilibrating groups of forces as the statically indeterminate 


A worked-out example shows how nine suiiably 


quantities. 
chosen groups suffice to caltulate the stresses in a wing that has 
7 statically indeterminate quantities according to the shear 
field theory. The method suggested was applied to the calcula 


tion of the stresses in a full-scale wing that was tested. Good 


agreement was found between theory and experiment 


INTRODUCTION 


— MEASURED AT THE ROOT of sweptback wings 
differ considerably from those calculated by means 
of the elementary theory of bending and torsion. The 
differences can be designated as the disturbance stresses. 
They are of great importance in the calculation of the 
ultimate load of the wing, as well as in the determination 
of its aeroelastic properties. 

When the shear field theory'~* is the 
stress analysis of an aircraft monocoque structure 


used in 
having m stringers and n ribs, a statically indetermi- 
nate problem of the order n(m — 3) must be solved. 
Usually m is large, and thus the calculation of the 
stresses is a highly indeterminate problem. The analy- 
sis is laborious, and, to reduce the work of computa- 
tion, it is general practice to replace the many actual 
stringers by a small number of imaginary longitudinals 
whose total cross-sectional area is equal to the effective 
bending area of the wing. 

When the stresses vary rapidly, this simplification 
often leads to an extremely rough approximation 
This is particularly true when pronounced stress peaks 
occur, and in such cases the stress values calculated 
may be too low. The method here described is pro- 
posed in order to reduce the computational work with 
out introducing such errors. It was developed for the 
stress analysis of sweptback wings, but it is also appli 
cable to cylindrical monocoque airplane structures. 


THEORETICAL CONSIDERATIONS 


The redundant wing structure is transformed into a 


statically determinate one if a sufficient number of 
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imaginary cuts are made in it in suitable locations. 
The necessary number of cuts is, as a rule, m — 3 at each 
rib station if the number of stringers is m. The forces 
caused by the external loads in this statically deter- 
minate fundamental system can be easily calculated. 
They are often called the zero forces. The actual forces 
in the complete redundant system naturally differ from 
the zero forces. They can be obtained if forces that are 
caused by groups of self-equilibrating forces acting in 
the cuts are superimposed upon the zero forces. It is 
known that the validity of the calculations is not im- 
paired if forces corresponding to the Bernoulli-Euler 
theory of bending and the Bredt theory of torsion are 
used as zero forces. This procedure is followed in the 
present paper. 

It is shown in the references* * how m — 3 elemen- 
tary self-equilibrating groups of forces can be chosen 
at each rib station as the redundant quantities and how 
another m — 8 more suitable groups can be obtained 
from them by means of linear transformations. These 
groups must be linearly independent. The criterion 
for linear independence is that the determinant of order 
m — 3 of the coefficients of the forces of the groups must 
not vanish 

When the structure is loaded with the forces of one 
of the m — 3 groups at a particular rib station, the two 
cells adjacent to the rib undergo deformations accord- 
ing to a well-defined pattern. An increase or a de- 
crease in the magnitude of the forces of the self-equili- 
brating group causes an increase and a decrease, respec- 
tively, in the magnitude of the deformations, but the 
deformation pattern remains unchanged. Each group 
of self-equilibrating forces can be considered a general- 
ized force, and the deformation caused by it can be 
denoted a generalized coordinate or a degree of freedom 
of motion of the system.‘ The multiplier of the uth 
group of forces is X,. It is calculated from the cond1- 
tions that the deformations of all the cells must be con 
sistent under the simultaneous action of the zero forces 
plus the n(m — 3) self-equilibrating groups, each 
multiplied by its proper factor X,. 
the gaps that open up between the various cells during 


In other words, 


the deformations caused by any of the individual groups 
of self-equilibrating forces or by the zero forces must 
close again when all of them, mixed in the proper ratios, 
act simultaneously. 

It is often possible to select a small number of degrees 
motion in such a manner that the re 


of freedom of 
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maining degrees of freedom of motion contribute only 
insignificantly to the actual deformations and internal 
forces. Then, the multiplying factors Y of these re- 
maining degrees of freedom of motion are negligibly 
small. They can be omitted from the calculations, 
with the result that the number of effective redundancies 
and the work of computation are greatly reduced. The 
process described is analogous to the Rayleigh-Ritz 
method. Whenever in doubt, the permissibility of the 
procedure can always be checked by considering a few 
more self-equilibrating groups in the calculation. If a 
pronounced stress peak appears in a longeron, it is ad- 
visable to increase the number of degrees of freedom in 
such a way as to load this part of the structure with 
particularly severe deformations (or forces). In this 
case, the calculated stress peak may be a little higher 
than that obtained from a rigorous calculation (this is 
in agreement with the Fourier theory), but the error 
cannot result in an unsafe structure. 

In a sweptback wing the cover plates are trapezoidal. 
The writer has shown earlier? how equilibrium can be 
established in such structural elements. 


NUMERICAL EXAMPLE 


The sweptback wing of Fig. 1, loaded in the manner 
shown in the figure, has been studied by means of the 
theory just discussed. In this wing, m = 12 and the 
number of redundancies is 9 X 3 = 27 when the load- 
ing is symmetric. Fig. 2 shows the stress distribution 
corresponding to nine generalized forces. The groups 
of forces are obtained by multiplying the mean stress 
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FIGURE 2 





at each longeron by the effective area. As the stres 
distributions are approximate, it is, of course, necessary 
to correct for the conditions of equilibrium. The nine 
generalized forces correspond to nine degrees of free- 
dom. Let us indicate the cell between the ribs i — | 
and i by 7. The zero forces are calculated according 
to the theories of Navier and Bredt. In cells V and 
VI, the zero forces must satisfy the following condi 
tions: (1) The forces in the longerons must have the 
same distribution (in the rib direction) as the Navier 
forces at rib station IV; (2) each panel in each cel 
must be in equilibrium. It is advantageous to choose 
the Navier-Bredt forces as the zero forces, because in 
this manner the disturbance stresses are obtained 
directly from the self-equilibrating forces after their 
X,, values are calculated. 

If a straight box beam is loaded with shear forces as 
shown in Fig. 2a, the corresponding stress distribution 
is approximately that shown in the figure.* If 4 
straight box beam with four longerons (Fig. 2b) 3s 
loaded by a torque, the stress distribution near the 
fixed part of the beam is approximately the one shown 
in Fig. 2b. When the box beam has more than four 
longerons, the resulting stress distribution is approxi 
mately the superposition of the distributions in Figs. 
2b and 2c. As a sweptback wing loaded by a sheaf 
force is subjected to bending, as well as to torsion, in the 

* Calculations have been carried out with several shapes 
(parabola and sine curves) and yielded small differences in the 
resulting stress distribution. 
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neighborhood of the wing root, the self-equilibrating 
groups of forces corresponding to uw = 1, 2, and 3 make 
the principal contribution to the disturbance stresses 
(Figs. 2a, 2b, and 2c). Fig. 3 presents the resultant 
forces at the critical section V of the sweptback wing of 
Fig. 1. 
result of a calculation with nh = 
In this case, the groups corresponding to 


The zero forces are given in curve a, and the 
1, 2, and 3 is shown 


in curve b. 
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SWEPTBACK WINGS 63 
uw = |, 2, and 3} were applied to rib stations IV, V, and 
VI, so that nine redundancies were used in the calcula- 
tions. The corresponding stress distributions are given 
in Fig. 4, curve b. 

For the calculation of curve c, the redundancies were 
increased by adding uw = 4 and 5 at rib stations V and 
VI. As the disturbance stresses are small at rib sta- 
tion IV, it was considered unimportant to apply u = 4 
and 5 to this rib station also. This calculation was thus 
carried out with 13 redundancies. 

In accordance with the stress function in Fig. 4, curve 
b, stress peaks exist at the flanges in the rear spar. It 
is important to know whether the stress values in these 
flanges are determined accurately enough with the nine 
redundancies listed. For this reason, another calcu- 
lation was carried out with n» = 6, 7, 8, and 9 (Fig. 2) 
applied at stations V and VI, together with u = 
and 3 at stations IV, V, and VI. The result is shown in 
As anticipated, the forces in the 


ve + 


Fig. 3, curve d. 
flanges at the rear spar are a little higher than the true 
values. The generalized forces corresponding to u = 
6, 7, 8, and 9 differ greatly from the true forces, and it 
is obvious that the true stress distribution cannot be 
obtained with these degrees of freedom. The calcula- 
tion shows, however, that the stress peaks in the rear 
spar are estimated with sufficient accuracy when the 
calculations are carried out with the nine redundancies 
chosen. 

Finally, Fig. 3, curve e, contains the result obtained 
LZ, 3, 4, 3, 6, ¢,.3, 
1, 2, and 3 


from a calculation made with up = 
and 9 applied to stations V and VI and yw = 
The corresponding stresses are 
were 


applied to station IV. 
shown in Fig. 4, curve e. 
carried out with 21 redundancies and prove, as can be 
seen from curves b and e in Fig. 4, that the use of nine 
redundancies suffices for all practical purposes. 

A calculation with the zero forces changed in pro- 
portion to the forces given in Fig. 3, curve a, was carried 
out (variations up to 12 per cent) with uw = 1, 2, and 3 
applied to stations IV, V, and VI. The curve obtained 
differs so little from curve b in Fig. 4 that the difference 
cannot be shown in a figure of this size. 


These calculations 
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Calculation of the Thermal Boundary Layer of 
a Body in Incompressible Laminar Flow 


W. Dienemann 

Dipl!.-Ing., Institute of Fluid Mechanics, Technical University, 
Braunschweig, Germany 

September 15, 1950 


Is A RECENT PAPER BY Goland,' the heat transfer on a certain 

infinitely long cylinder has been calculated exactly and com- 
pared with the results of several approximate calculations, which 
gives good agreement only in one case. The exact calculation 
is based on a solution by Sears? for the velocity distribution in the 
spanwise direction of the boundary layer of an infinitely long 
yawed cylinder, which leads to the same differential equation as 
the thermal boundary layer if Prandtl Number Pr = c,u/k = 
v/a = 1, (a = k/pcp). On the same lines, an approximate solu- 
tion for the thermal boundary layer might be given according to 
Wild,* who has applied the Pohlhausen* method to the velocity 
boundary layer of the infinitely long yawed cylinder. This has 
already been tried by Kroujiline® and Fréssling.® 

We want to explain here that this approximate solution for the 
thermal boundary layer can be simplified considerably if the 
Holstein-Bohlen’ version of the Pohlhausen method is applied. 


On integration of the energy equation 


u(OT/dx) + v(OT/dy) = a(d?7/dy?) 


with respect to y (coordinate normal to the surface of the body), 
one gets the following ‘“‘heat-flow equation’’ of the thermal 
boundary layer: 


7 


d ; od, * 
(ud,*)\dy = -a( ) 
dx 0 oy 0 


where 3;* = (T — 7To)/(T. 
ture distribution in the boundary layer. 
the edge of the boundary layer; 7, = temperature of the wall.) 
The heat-flow equation |Eq. (1)] corresponds to the momentum 
equation of the velocity boundary layer. 

We now, in analogy to Holstein, introduce the ‘‘heat-loss- 


(1) 


Ty) is the dimensionless tempera- 
(7) = temperature at 


thickness” of the thermal boundary layer. 


u 
v7: = x i dy 
J 0 l 


is the velocity distribution in the boundary layer 


(2) 


where u/U 
Putting 


6= (3) 


Vet? /v 


one gets from Eq. (1) the following ordinary differential equa- 


tion for 6(x): 
(4) 


do/dx = F(A, A; Pr)/U 


which can be solved by the isocline method. [In Eq. (4) A; 
means thickness of thermal boundary layer 6,/thickness of ve 
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locity boundary layer 6,; A = (6,2/v)(dU/dx) = parameter alter 
Pohlhausen.] The main advantage of our method js th 
F, (A, &, Pr) isa universal function not dependent on the Speci 
problem. By the method of Wild (compare also Eckert), 9; 
would get a much more complicated differential equation 4, 
Ai(x), with coefficients different for each special case, the sSoluti 
of which can be arrived at only by laborious calculations 

Integrating Eq. (4) yields 6(x) along the surface of the body 
furthermore, by Eq. (2), the thickness of the thermal boundap 
layer &, because u/U is known from the Pohlhausen method ap 
the temperature profiles 3;* are assumed to be affine and give 
by acertain polynomial in y/&, 


’ ,\3 , 4 
i*=1-27 +2 ”) - ”) 
FY 5 5, 


In this way, also, all other boundary-layer parameters can } 
calculated and, especially, the rate of heat transfer (local Nussel 
Number). 

Our approximate method works satisfactorily for the therm 
boundary layer with constant wall temperature and moderat 
velocity (generation of heat by friction neglected) but with arh 
The r 


sults of our method are in excellent agreement with experiment 


trary shape of body and arbitrary Prandtl Number. 


and with exact solutions (compare, e.g., Goland, reference 
Fig. 1). Eckert’ has given another approximate solution for th 
same problem which also yields good results and which is st 


simpler to carry out. 
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Fig.1: Heat transfer on a cylinder. 


On the other hand, our method has the valuable advantage that 
it can be extended to other more general problems of therma 
Up to now this has been done for two cases 
|taking into account the heat 


boundary layers. 
(1) the “thermometer problem” 
generated by friction and assuming an insulated surface, no heat 
transfer, (O7'/Oy)y=0 = 0}; (2) the heat-transfer problem for 
cylinder with arbitrary variable wall temperature 7\(x) 

For the thermometer problem, there are no important altera 
tions of our method as compared with that given above. For 
the heat-transfer problem with variable wall temperature, the 
temperature profiles in the boundary layer can no longer be as 
sumed to be affine. Eq. (5) for the temperature distribution has 
to be extended in the following way: 


; T — Tp é (2) (?) 
v1 = a m I(x)? FT - + Ar: G (O 
AT» 5: 5: 


where F and G are universal functions. Here, the wall tempera 


ture is given by 


T w(x) — To = AT: f(x) 
AT) being constant, and the new parameter A; is proportiona 
to df/dx = f'(x) and given by 
Ae = (6:2/v) Ap: Pr- U[2 + (A/6)]- f’(x) 7 
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In this case one gets, for 6(x), the differential equation 


d6— Fa f(x), Ae A, Me, Pr] 


dx U 


(8) 


where again F, is a universal function of the parameters indicated, 
not depending on the special problem. Eq. (8) again can be 
solved by the isocline method, thus giving all parameters that 
are needed. 

As an example of Eq. (8), the thermal boundary layer has been 
calculated for 2 flat plate of length L, heated at the front half, 
and cooled at the rear half with f(x) = cos (rx/L). In agree- 
ment with Chapman and Rubesin,’ we got the result, at the first 
moment very surprising, that at that place where the temperature 
of the plate is equal to that of the flow outside of the boundary 
layer, x = L/2, there is heat flowing into the plate. The solution 
of Chapman and Rubesin is restricted to the flat plate and 
needs an approximation of the wall temperature 7,,(x) by a poly 
nomial in x, but it takes into account compressibility and, there 
fore, is extremely complicated. On the contrary, our method, as 
outlined above (neglecting heat due to friction), is applicable to 
an arbitrary shape of body and to an arbitrary steady distribu 
tion of the wall temperature but is, nevertheless, simple. It is 
expected that our method can be extended still further and that 
it will be possible to examine separately all different parameters 
that have an influence on the thermal boundary layer. 

jur method, as outlined here briefly, will soon be published in 
detail as a doctoral thesis of the Technical University of Bruns 


wick, Institute of Fluid Mechanics. 
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NoteTon ‘‘Longitudinal Stability of Autopilot- 
Controlled Aircraft’’ 


G. J. Sissingh 
Royal Aircraft Establishment, Farnborough, Hants., England 
September 20, 1950 


I THE PAPER ‘“‘Longitudinal Stability of Autopilot-Controlled 
Aircraft’? by A. Vazsonyi (JOURNAL OF THE AERONAUTICAL 
SCIENCES, July, 1950), it has been assumed that the signal and 
the elevator motion are connected by a linear ordinary differen- 
tial equation of the first order—see Eq. (94) of the original paper, 
which, using the same notation, can be written as 


6+ 7,5 = 70 + uT46 (1) 


FORUM 65 


Similar simplifications have been made previously by other 
authors. It may be of mention that mechanical 
control devices of this kind do exist and that the Young-Bell 
stabilizer of the Bell helicopter and the servo blades of the 
In these cases, however, 


interest to 


“Hiller 360”’ work along these lines. 
This means that the servo signal has only a derivative 


z= 0. 
term. 

I should like to add another remark that gives some further 
information on the physical behavior of control devices char 
acterized by Eq. (1). As can be seen from the representation of 
stability problems in the vector diagram, it is not the time lag 7, 
(i.e., the time expressed in seconds) which really matters but the 
phase of the control application or, more accurately, its phase 
The phase angle of autopilots with time lag as character 
(1) depends on both the fre 


angle. 
ized by the left-hand side of Eq 
quency and the damping of the oscillatory motion of the aircraft 
The latter effect, however, is of minor importance and can there 
fore be neglected. For a pitching oscillation with constant amp 
litude the phase angle ¢ is given by 


tang = 7\w (2) 


where w is the frequency of the oscillation in pitch. This means 
that the time lag 7, has the effect of an apparent change of the 
control characteristics. Or, in other words, for a given fre 
quency w, the autopilot with the time lag 7); applies the same 
control displacements as an ideal autopilot without time lag 
having the control characteristics 

P — a+ ww? TT, ul, — pT, (3) 


1+ wT 1+ 7) 


By comparison with Eq. (1), it follows that the time lag 7) 


changes the proportional term from { to 


(f + pw?T47,)/(1 + w?T;?) 


and the derivative term from a7‘, to 

(ul 4 = aT) (1 + w*T,?) 
If several natural frequencies occur, Eq. (3) can be applied to 
each of them. This means that for any mode of oscillation the 
autopilot can be replaced by an individual ideal control device 
without time lag that results in the same oscillation of the air 
craft. 

In complex (preferably neutral, slightly stable or unstable) 
systems with control devices, the input signal of which cannot 
be arbitrarily changed, this method is sometimes useful to find 
out which modifications of the control device or the system are 
necessary for improving stability. This is done by (a) comput- 
ing the control characteristics of an ideal autopilot without time 
lag, which results in a satisfactory stabilization of the aircraft, and 
(b) comparison of the desirable ideal control characteristics with 
those obtained by Eq. (3) with w = frequency of that mode of 


oscillation which is to be stabilized 


On the Solution of Total Differential, 
Boundary Value Problems 


Arthur N. Tifford 

Associate Professor, Aeronautical Engineering Department, 
The Ohio State University, Columbus 

September 28, 1950 


N 1943, WHILE PERFORMING the numerical integrations that 
I were reported in reference 1, the author employed a simple, 
general procedure that enables the solution of two-point bound- 
ary value, linear, total differential equations to be obtained di 
rectly without recourse to trial and error. Because of its sim 


plicity, he assumed the procedure to be well known. During the 
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last year or so, however, he has discovered that it is not. Ac- 
cordingly, this note has been written to bring it to the attention 
of all. Some remarks as to its significance in the solution of non- 
linear differential equations have been added. 

The procedure is most easily understood by following its ap- 
plication to a relatively simple case. Consider the second- 
order differential equation 


VY" + fi(x)V’ + fo(x) Y = fa(x) 
(atx =a, Y = 5b; atx =c, Y=d) (1) 


where the functions, fi, fe, and f;, are given, known functions 
and a, b, c, and d are known values. It is desired to find the 
value of Y’ atx = a [writtenas Y’(a)]. 

We first obtain the value of Y and its derivatives at the first 
few steps taken, starting at one boundary toward the other. 
This may be done by expanding the unknown function, ¥, as a 
Taylor series near the starting boundary. 


: ; a (Az). Ci 
Y = Y(a) + AxY’(a) + a Y’'(a) + = Y’’’(a) + 
(Ax)* __ 
Y''""(a) +... +, ete 
4! 
(2) 
" i See (Ax) | ,,,, 
Y’ = V’'(a) + (Ax) Y"(a) + ~ Y’"’(a) + 
Cl or 
31 } (a) +... +, ete. 


The values of the ¥"(a)’s are found as follows: One boundary 
condition states that V(a) is b. Y’(a) is unknown. From Eq. 
(1) we obtain ¥’’(a) as 


f(a) — fila) Y'(a) — feo(a) V(a) 
or 
Y’’(a) = Fo(a) + Y’(a) Fr(a) 


The value of ¥’’’(a) is obtained from Eq. (1) after the latter has 
been differentiated. Thus, 


Y’"'(a) = F;(a) + Y’(a) F33(a) 


Similarly, the higher derivatives are found by successive dif- 
ferentiations of Eq. (1) as 


Y"(a) = Fi(a) + Y'(a) Fan(a) 
Substitution back into the Eqs. (2) yields 
Y™(x) = F(x) + Y’(a)Fmm(x) (3) 


for Y and its derivatives for the first few steps of integration. 

We are now ready to apply any one of the available methods 
of numerical integration (the author prefers the Milne method?) 
to a step-by-step approach to the outer boundary. Once there, 
we require Eq. (4) to be satisfied. * 


d = Fi(c) + Y'(a) Foc) (4a) 
or 
Y'(a) = [d — Fy(c)|/Foo(c) (4b) 


The procedure just described is the one used by the author 
many times in manual numerical integrations of boundary-layer 
equations. It is, of course, directly applicable to digital machine 
computations. It is not useful, however, in connection with 
analog computer determinations. For the latter, the following 
related approach has been developed: 

Suppose that, instead of carrying ¥’(a) along as an unknown 
in our analysis, we assign an arbitrary value, A, to it. Then, at 
the second boundary, we have Eq. (5) instead of Eq. (4a). 


Ya(c) = Fo(c) + A Foo(c) (5) 


* In a boundary-layer type of protlem, c is infinite and therefore unat 
tainable. In such a case, however, the value of Y at large values of x is 
essentially the same as at infinity. We therefore apply Eq. (4b) (in which 
x has been substituted for (c) at several successive, large values of x When 
the values of Y’(a) so determined do not change significantly with x, we 
have the value of Y’(a) sought 


where Ya(c) differs from the required boundary value, d, 
Fy(c) and Foo(c) are unknown. 
Now, assign a second arbitrary value, B, to Y’(a) and integray, 


Eq. (6) is found applicable at the second boundary. 
Ya(c) = Fo(c) + BFool(c) 


where Ypa(c) differs from the required boundary value, d, a 
Fy(c) and Foo(c) are the unknowns of Eq. (5). 

Ya(c), Ya(c), A, and B being definite numbers, Eqs. (5) a 
(6) constitute a set of simultaneous equations from which th 
values of Fo(c) and Fyo(c) can be determined.t These values, 
Fo(c) and Foo(c) are then substituted into Eq. (4b), and th 
proper value for }’(a) is determined. 

The advantage of these methods is more fully apprecia 
when applied to more complex problems. For example, jn 
current application, two simultaneous fourth-order equations ; 
two dependent variables with four boundary values specific 
at each boundary are handled. Because of the four simultane, 
unknown conditions at each boundary, solution by means of 
cut-and-try procedure is a herculean task—even when using y 
analog computer. The methods described, on the other hap 
simplify the solution procedure to essentially five integrations 
done simultaneously, if desired—and the subsequent solution ; 
four simultaneous linear algebraic equations. 

It is stimulating to consider the application of the same ide 
to the solution of nonlinear total differential equations. A litt 
thought on the matter reveals that, in the strictest sense, th 
above linear algebraic equations at the outer boundary are r 
placed by equations of infinite series. (If the step-by-st 
integration were actually carried out, however, we would obtai 
equations of finite, high-powered, polynomials instead becaus 
of the approximation involved in using steps of finite size 
Thus, it is apparent that there are a large number of mathemati 
solutions to a given nonlinear problem. Of course, in genera 
most of these solutions are complex or imaginary, but there 
no guarantee that only one real solution exists. In other word 
whereas we can be sure that the solution obtained to a linear dif 
ferential equation will agree—barring numerical errors—wit 
the physical solution, in the case of nonlinear equations, add 
tional considerations (such as the stability of different solution: 
must be brought to bear, in general, before we can be sure of suc 
agreement. 

Because of the relative directness of the numerical solution 
linear differential equations, it would appear advantageous t 
solve complex nonlinear differential equations by a procedur 
involving linearization followed by successive perturbation 
The relative merits of this and related procedures are to lk 
studied shortly. 


+t In more complex problems having many unknown initial boundar 
conditions, a judicious choice of values for A, B etc., greatly simplif 


the set of simultaneous equations 
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Concluding Remarks on the Speed of Sound 


M. Morduchow and M. Morkovin 

Polytechnic Institute of Brooklyn and University of Michigan, 
Respectively 

September 27, 1950 


gears 1 AND 3 POINTED OUT the mathematical diff 
culties in the customary textbook derivations of the for 
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for the speed of sound 


mula 


(dp/dp) (1) 


9 
eG = 


based on the steady one-dimensional equations of nonviscous 
quids. The difficulties center around the interpretations of 
differentials du, dp, dp, etc., occurring in the momentum, energy, 


ind continuity equations. The usual interpretation employed 


in the derivations of these equations of motion—namely, 
du dp dp 
du = dx, dp = dx, dp = dx 
dx dx dx 


with the limiting process dx — ( in mind does not appear 
well older 


etc., 
rigorously applicable. 
references such as reference 4) have indicated that these difficul- 


References | and 2 (as as 
ties disappear when the sound wave is considered as the limit of 
a shock wave as its strength approaches zero. Thus, if a shock- 
strength parameter is considered as the independent variable in 
the limiting process (rather than x) and the increments across 
the shock waves Au, Ap, etc., are related to this parameter, the 
algebraic manipulation of the usual derivations become mean- 
ingful and rigorous. (It is believed that this limiting process 
is implied, though without explanation, in some of the textbooks, 
such as reference 5. ) 

One purpose of this note is to present clearly the details of 
this approach which are not usually available to the novice in the 
field of compressible flows. A second purpose is to present to 
the same public (which may have been somewhat coafused by 
some of the disagreement between references 1, 2, 3, and, more 
recently, 6) the wide area of agreement reached through these 
discussions. It is felt that the Readers’ Forum will further in- 
crease its present high usefulness to the general public if the 
positive results of any discussion or controversy are stated and 
emphasized. It is hoped that a precedent may thus be estab 
lished. 

Consider a discontinuity (stationary with respect to an ob 
server) separating downstream flow characterized by wt, f», 
p,, and 7», from upstream flow characterized by ™, pi, p1, and 
T;. The fact that the discontinuity is stationary implies that 
its speed of upstream propagation with respect to the fluid, a, 
is numerically equal to m. The feasibility of this one-dimen 
sional flow system rests completely on the consistency of the con- 


servation laws and the equation of state. 


pill, = pele (‘> 

pyuy? + hi = p22” + pe (M) 
(1/2)uy? + CpT, = (1/2)u2? + ClT2 (E) 
pi = Rol; po = RooT> (S) 


Setting vw. = u, + Au, po = pi + Ap, p2 = pi + Ap, and elimi 
nating T from Eq. (E) by the use of Eq. (S), the conservation 
laws can be written in the useful nondimensional incremental 


form 
Au Ap ] ” 
+ = () i") 
My pi 1 + (Ap/p1) 

(Au/u) + (pi pity?) (Ap pi) = 0 (M’‘) 

Ap At Au \? 2 A A 
fi+ ‘ “+ ( ‘ 7 n) Pp _ 4e| _, 

a pi uy uy y — 1 \pim,? pi pi 
(E’) 


The conservation laws thus form a system of three equations in 
the four nondimensional quantities: p;/p:m?, Au/m, Ap/p,, 
and Ap/p;, all four of which are measures of the strength of the 
discontinuity. The system of equations is consistent, validating 
the original assumptions. Furthermore, one can arbitrarily 
assign values to one of these quantities (i.e., choose it as the 
strength parameter), say (Ap/p,), and expect that the other three 
will be determined from Eqs. (C’), (M’), and (E’). This is 
indeed the case; (Au/u,) is immediately determined from Eq. 
(C’). Substituting for (Au/u,) and (Ap/p,) from Eqs. (C’) and 
(M’) into Eq. (E’) yields, for Ap/p; # 0, 
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Ap — | 
e=]i- x of & (2) 
Pl pili,” 2 py” 


Eq. (2) determines the quantity (f)/p,“,?) for every chosen 
value of Ap/p;. Finally, with (p;/pyu:?) and (Au/u,) known, 
Eq. (M’) determines (Ap/p,). Thus, as the strength parameter 
Ap/p; varies, the nondimensional characteristics of the flow sys- 
tem under consideration are completely and uniquely deter- 
mined. Mathematically, a sound wave can be defined as a limit 
of shock waves for which the strength parameter approaches 
zero, Ap/p, > 0. This limiting process applied to Eq. (2) re 
sults in the determination of the speed of propagation of sound 


ui? = a? = y(pi/p1) (3) 


(3) leads to a more familiar interpretation of the non 
dimensional quantity p,/p.m:7—namely, 1/yM,? M, is 
the upstream Mach Number. Incidentally, Eq. (2) now pro- 
vides a useful formula for the compression Ap/p; in a finite shock 


Eq 
where 


wave 
Ap/p, = (M2 — 1)/f{1 + [(y — 1)/2)43} 
Subtracting Eq. (C’) from (M’), one obtains 
ui? = (Ap/Ap) [1 + (Ap/p:)] 


Considering, as before, the limit of a weak shock wave and thus 


letting Ap/p,; —~ 0, it is seen that for a sound wave u,? = 
lim (Ap/Ap). If, as explained previously, the quantities po, 7», 
Ap —> 0 


and uw in Eq. (C), (M), (E), and (S) are considered as functions 
of a shock-strength parameter Ap/p:, then the limit of the ratio 
can be interpreted as an ordinary derivative. Hence, 


(dp/dp) (4) 


uy? = 


To show, by the weak-shock-wave approach, that a sound 


wave is isentropic, the entropy change, 


Se — S: = Cy log (2 (":) 
Pi/\p2 


across the discontinuity can be calculated in terms of the strength 
parameter Ap/p;. Elimination of Au/m, and p;/p,u,? between 
Eq. (C’), (M’), and (E’) leads to the well-known Hugoniot rela- 


/ 
pe y¥ +1 Ap y — 1 Ap 
= 1+ = / i— * 
pi 2 pi J | “ pi 
Since AS in Eq. (4) can now be expressed as a sum of terms of the 
form log(1 + const. Ap/p,), it can be easily expanded into a series 
that starts with the cubic term 


(5) 


AS = 


tion 


(6) 


: C, ; Ap \* Ap \' a 
AS = —-y(y? — 1) + 0 (7) 
12 pi pi) 
Hence, 
lim (AS/Ap) = 0 
Ap > 0 


and a sound wave viewed as the limit of a weak shock is isen 
tropic. 

It is thus seen that the familiar expressions (3) and (4) for the 
velocity of sound, as well as isentropy, can be derived by a simple 
means without any of the apparent mathematical inconsistencies 
which a direct differential-equation approach would contain. 
It should be noted, moreover, that the approach indicated here 
namely, a shock wave becoming weak to 
rather clearer physical picture than the situation implied by the 
use of the differential equations. 

It is clear that a full presentation of the preceding lengthy an 
alysis in a classroom may, if desired, be omitted if the essence of 


seems present a 


the argument is outlined. 

Professor Clark B. Millikan has indicated to us that the pre 
ceding treatment is conceptually equivalent to the one he uses 
and that it is implied in essence in his recent note.® 
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On the Calculation of the Potential Flow 
Around Airfoils in Cascade 


N. Scholz 
Dr.-Ing., Institute of Fluid Mechanics, Technical University, 
Braunschweig, Germany 


September 15, 1950 


N ORDER TO IMPROVE the aerodynamic properties of airfoils in 

cascade, it is necessary to have a convenient method for the 
calculation of the potential velocity field. The method of con- 
formal -transformation, in-our opinion, is much too complicated 
for practical purposes; only the ‘‘singularity method,’’ where 
the airfoil is replaced by a certain distribution of vortices and 
sources and sinks, seems adequate. In some recent papers,'~* 
the singularity method has been applied to airfoils in cascade. 
The calculations on these lines can still be simplified considerably, 
as we have shown in a paper’ to be published shortly in detail 
and an extract of which may be given here. Our method has 
already been applied with good success to the calculation of 
boundary layers on airfoils in cascade.* 

According to the singularity method, each airfoil of the cascade 
is replaced by a continuous distribution of vortices and sources 
(and sinks), situated on the mean line of each profile. The main 
problem is to calculate the field of induced velocities due to these 
singularities. 

In the complex s = x + /y-plane (Fig. 1), an infinite row of 
elements of singularity of strength g-dz’ may be given situated 
on a straight line through the origin inclined at the angle \ 
against the y-axis. The distance of the airfoils being ¢, these 
elements have the complex coordinate 


zy =twte™ (vy =0, +1, +2,...) 
In the general case, the strength g is complex, g(z) = g(s) + 


iy(s), g(z) meaning the distribution of sources (and sinks) and 
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Fig 1: A cascade of airfoils in the complex plane z= x+iy. 


y(z) meaning the vorticity. The field of induced velocities du 
toa single element g(z’)- dz’, situated at 2,, is 
d 


dw,(z) = du, — idv, = g dz'/2n(z — wte~*>) 


Summing the contributions of all elements of the infinite row 
yields 
, +@ 
g dz 1 
dw(s) = =a = 
Zte~** = —ao (wz/te~**) — tay 
This infinite sum is identical with the series expansion of hype, 


bolic cotangent with the argument e’xz/t, so that Eq. (1) ex 
] Ca) 


be written 


ets , TZ Yi 
dw(z) = g(s’) coth e'“ } dz la 
2t t 
Summing up overall elements of the mean line, 2’ = ~¢/2; 


+ c/2, finally gives the induced velocity field of the whole cascay 


etd +e/2 ee, 
u(s) = g(s’) coth { we™ Jae 9 
2t 2’=-¢/2 t 


Here, \ is identical with the angle between the direction of th 
mean flow U..* along the x-axis and the normal to the cascad 
axis (Fig. 1). 

For airfoils with small camber (s’ ~ x’), and taking the ip 
duced velocity at the chord only (x-axis), one gets from Eq. (2 


46/2 
eth 


i p eo id , 
w(x) = g(x’) coth { ze’ dx 
2t x'=-—c/2 t 


in analogy to the downwash on an isolated airfoil according 
the theory of Birnbaum and Glauert.' In the limit t > ©, Ry 
(3) simply gives 


W(X )isol. = - ax" 38 


the well-known formula of an isolated airfoil. In the lim 
t— 0, Eq. (3) can be put into the form 


r +c/2 ’ 
g(x’) 
2t J x’ c/f2 t 


where I is the tctal circulation of one airfoil and I°/t is the tots 
circulation per unit length of the cascade axis. 
For the numerical evaiuation of Eq. (3), the singularity of t! 


integrand at x = x’ is inconvenient. This singularity is a matte 


of the isolated airfoil only, and therefore it is reasonable to spl 
off the induced velocity field of the isolated airfoil from the ve 
locity field of the cascade. Thus, one gets from Eq. (3) 


W(X )easc. —°W(X )isol. = 


, . i= ¥" t 7 
g(x’) | e'* coth [ re" = 1d 
2t x’ c/2 t w(x — x 


The numerical evaluation of Eq. (4) is convenient because th 


bracket |] is regular in the whole range of integration and is 


universal function of (x — x’)/t and \ only, not dependent on th 


special airfoil. This universal function has been tabulated i 


reference 7. The evaluation of the velocity field of the isolate 


> 


airfoil, Eq. (38a), can be done very simply by the method « 


Allen. In this way, for a prescribed distribution of vorticity 


and sources, the induced velocity along the chord of the cascaé 
airfoils can be easily calculated 


From Wease. = Uease. — (ease. and with the superposed trams 
verse velocity U. along the x-axis, one gets the total velocity 


the chord c of the cascade airfoil 


W(x) = V (Ua + theaac.)* + Veasc.” 


* This is the geometric mean value of the directions of flow far in fros 


and behind the cascade. 
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The shape of the profile y(x) is obtained from the kinematic flow 


condition 
— Kiso.) | (6) 


dy/dx = Ucasc. [U + (Ucase 


Finally, the velocity distribution on the contour of the profile 
Pis obtained from that on the chord according to Riegels® by 


We(x) = Welx)/W1 + (dy/dx)? (7) 


which yields a stagnation point (Wp = 0) at the leading edge. 
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A Method of Correction of Astigmatism in 
Schlieren Systems 


Rochelle Prescott and E. L. Gayhart 

Applied Physics Laboratory, The Johns Hopkins University, Silver 
Spring, Md.* 

October 10, 1950 


I CONNECTION WITH some experimental work with the schlieren! 


method for the observation of airflow, a simple economical 


method of correcting the astigmatism in the common two 


mirror system by means of a planocylindrical lens has been found 


In this work it was desirable to use a ‘“‘point”’ source and a “point” 


stop or knife-edge in a system that consisted of the common 


Z arrangement of two ordinary paraboloids as shown in Fig. la 
The astigmatism was of the magnitude indicated in Fig. le. As 


a result, it was possible to use only a vertical or a horizontal 


knife-edge in the position of the corresponding stigmatic image 
The difficulty was resolved by the use of a positive planocylindri 
cal lens of the proper power between the light source and first 


mirror as shown in Fig. la. To date, no deleterious effects have 


been noted 


The separation of the tangential and sagittal images was 


and the value of lens was computed which at a reason 


measured, 


able distance would give an image displaced by this amount. 


In this particular case, it turned out that a three-quarter diopter 
lens would give the required displacement at approximately 150 
mm. Such a lens was procured and mounted so that final ad- 


justments could be made by trial and error. The final adjust 


ments were made by using a microscope to observe the source 


image 
The results were remarkably good 
artificial star 


A critical examination of 
a West 
that the 
The pattern, however, 


the Fraunhofer pattern of an (an image of 


ern Union 2-watt concentrated are lamp) did show 


resolution of the system was not perfect 


was good. In actual use as a schlieren system with a Gaertner 


Spectrograph Slit set at 0.03 mm., the field could be darkened as 
uniformly after the introduction of the correcting lens as without 


* Operating under Contract 7386 with the Bureau of Ordnance, U.S 
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Fic. 1. Diagram of schlieren system showing relative positions: 
CORRECTED 
UNCORRECTED 
Fic. 2. Photographs of images at focal plane. (Magnification 
of image = 8 
it. It was possible to use either a vertical or horizontal knife- 
edge at the same position along the optical axis. It was also 


possible to use a ‘‘point’’ source as desired and have the field 


uniformly illuminated. 

In Fig. 2 are shown photographs of successive images at the 
focal distance of the second mirror. The positions of the images 
Ib and le. Comparison of the photographs 


after correction, 


are indicated in Figs 


with the diagram shows that, the critical image 


of the pinhole source is formed in the same plane as that in which 


the pinhole is imaged as a vertical line before correction 


As a result of the success of the above correction, the image of 


the field in the schlieren plane was examined critically with a view 


to correcting it for astigmatism by means of a second lens placed 


knife-edge and camera lens, if desirable rhe 


however, proved to be so small that it was felt cor 


rection was not worth while 
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An Approximate Treatment of Multishock 
Bodies of Revolution at Small Angles of 
Attack by the Method of Characteristics 


R. E. Bower 

Research Engineer, Grumman Aircraft Engineering Corporation, 
Bethpage, L.!., N.Y. 

October 17, 1950 


y BECOMES EVIDENT from an extended treatment of the bound- 
ary conditions at shock fronts that the method of character- 
istics as employed by Ferri! is applicable in certain instances to 
bodies of revolution generating more than one shock surface. 

The basic theory assumes that cylindrical velocity compo- 
nents u, v, win the coordinate directions x, y, ¢, respectively (Fig. 
1), are a solution of the fundamental relations governing the 
motion of a steady and inviscid supersonic medium when ex- 


panded as 


u = uy, + ale COS ¢Y 
v =v + atrcos ¢ ? (1) 
w= aw, sin ¢ 


where %, us, 1, ..., are functions of x and y alone and where 
the angle of attack of the body, a, is considered small enough to 
permit neglecting terms of order a?. This requires that Eqs. (1) 
satisfy boundary conditions at shock waves as established by 
Rankine and Hugoniot and necessitates the proof of the existence 
of shock surfaces deforming the flow field in a way commensurate 
with Eqs. (1). The present note is concerned with such a proof 
and is essentially a generalization of the proof given in reference 
1 for the existence of a single shock surface emmanating from the 
nose of a body. 

The shock surfaces are assumed to be described as follows: 
A body of revolution at zero angle of attack generates curved 
shock surfaces of revolution which can be considered to be the 
envelopes of right circular cones with axes coincident with the 
body axis and generatrices as lines tangent to the shock surfaces. 
For the body at an angle of attack, the resultant shock surfaces 
are no longer surfaces of revolution but are, for small incidence 
angles, assumed describable as the envelopes of the same cones 
considered when a = 0 with their axes rotated through small and 
variable angles, 7, measured in a ¢ = 0, w meridian plane. 

A representative “‘tangent shock cone”’ is shown in Fig. 1 for 
a =UVanda #0. The equation of the yawed cone is 
tan—! (y/x) = n+ cos ¢ (2) 


where » = constant = semivertex angle of the cone. A vector 


normal to this surface is grad.[tan~! (y/x) — uw — ncos ¢], witha 
unit vector written 
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Outward along a conical ray, the unit vector takes the form 
oe 
2 
§r: | ee 
Vixr + y? Vat + y? 
ae —s 
hence, the unit vector normal to the plane formed by &y and &7 is 


i. a oe 
tw =&rxX Ey or 
re . x . 
Ew: nsing, — ~ sing, 1 
y 
sti 


The desired breakdown of the velocity vector V at a point P re- 
sults in a component normal to the a # 0 shock cone as 
—_ 
Vy = —V-&y = usin» — veosu + 
ncos ¢ (vsinu + ucos yu) (3a) 


and tangential components 
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Vr = Vr = ucosut+vsing + 
n COs ¢ (Vv COS w — usin y) (3p) 
and 
Vw = V-éw = w+ sin 9 (u — vcot pz) 30 


Eqs. (3) are rearranged from a consideration of Eqs. (1), wher 
the quantities with subscript 1 are known from the a = 0) solu. 
tion at P(x, y). While these quantities are independent of » 
for given x and y coordinates, points on the a # 0 shock fron; 
[Eq. (2)] have variable coordinates x and y with changing ¢. |; 
is therefore necessary to express the components w and 9 as 
function of the flow properties at a point Q, near P, having fixed 
x and y values for all ¢. Specifically, Q is taken on the qa = (j 
tangent shock cone associated with the point P (see Fig. 1), j 


follows that 


On; 7 
Up = |\%m) Q + aN lo AN (4a 
and 
Ov; , 
MN\p = |%) Q + an % AN (4b 


where the derivatives are evaluated in the meridian plane ¢ = 
constant (plane OAQ) and AN is the distance QP’ given as 


AN = (|x! @/cos wu) n cos ¢ 4 


Defining the following quantities: 


Vy, = &) Sin wp — % COS wu 
Vr, = m cos uw + 9, sin uw 
Vivo = U2 SiN w — Ve COS wu 
Vireo = ue cos wp + v2 Sin wu 


substitution of Eqs. (4) and (1) into Eqs. (3) yields 


Vy = | Vy, . + acos¢ Valo + ncos g| Vm4|, + 


Q 
XiQ re) V7; 
n COS ¢ : (5a 
COS ON Q 
Ve = Viri\g + acos ¢|} T2\g — COS Vail + 
|X\Q oVr, , 
n COS ¢ = ob 
COS pu ON @ 
and 
Vw = asin ¢/we\g + 1 sin g(a, — % cot ue (5c 


> “qe . — . . 
The shock equilibrium conditions are introduced for axially 


symmetric flow as 





Vir,’ = Vr,” 6a 
Vu,"Vuy' = (Vmaz.2 — Vr,"?) (6k 
> ee a 
Vr 
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and for unsymmetrical flow phenomena 


Vr’ = Vr" (6c) 
Vw’ = Vw’ (6d) 
oe ee, ak eee - ' 
Vv" Vy = ( Vez.” V7"?) (6e) 
y+1 


where 7 is the ratio of specific heats, single prime indicates values 
in front of the shock wave, double prime indicates values behind 
the shock surface, and Vinar. is the maximum flow velocity real- 
ized when the fluid is allowed to expand adiabatically into a 
yacuum. 

Treating the flow entering or leaving the shock surface to be of 
the general form of Eqs. (5), substitution into Eq. (6c) yields, to- 
gether with Eq. (6a), 


‘ — - n \x0Vr," 
V1 = [V1 as J Ni a 
Q e a Q acosu| ON ¢@ 
“ee n |\x0Vr,"" 
| Ni — : (fa) 
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Ina similar manner, substitution of Eqs. (5), (7a), and (6b) into 


Eq. (6e) gives ‘ 
; Wins n  |xdVy," 
Vive - Vir = , = 
@ a e a@ COS u ON Q 
V2’ Vy,” . Ui Vr,' Vn," n xVy,”" OV y,’ 
Vy,’ 9 a Vy,’ a@ COS pw Vu,’ ON lg 
Ay—1)| Vr'Vn" 29 Vn" mxVr," OV 7, “1 
5" = T1 = ; : (4D) 
¥ + ] if vy,” a a COS ul V1 ON Q 
whereas substitution into Eq. (6d) shows that 
w"|g = | —(n/a) (u," — m4” cot w) + we’ + 
(n/a) (u' — %' cot w)\g (7c) 


The ratio y/a@ appearing in the above relations is independent 
of a for small angles of attack. To demonstrate that Eqs. (7) are 
also independent of ¢ and therefore show that the flow leaving the 
shock surface is in accordance with Eqs. (1), certain specific 
assumptions must be made regarding the type of flow entering 
the shock wave. 

For the initial shock wave off the extreme nose of the body 


(case treated by Ferri), the entering flow is constant and uni- 
Co 
formly inclined to the body axis. For a velocity Vy at an angle 


of attack a, 


u’ = Vo, v’ = aVocos¢ 

w’ = —a Vosing, Vy,’ = Vosinu 
Vy.’ = Vo cos pw, Vr,’ = Vocos u 
Vr.’ = Vosin yu, we’ = —Vo 


The ab- 


sence of y in Eqs. (7), after the proper substitutions, proves the 


with all derivatives with respect to N identically zero 


compatibility of the assumed shock surface and Eqs. (1). 

If a second shock-wave surface is present in the flow field, the 
entering conditions would, by reason of the above, be given by 
Eqs. (1), where no terms are zero but where #2, v2, and w» are inde- 


pendent of a and ¢. 
a shock surface remains unaltered in form. 


Eqs. (7) disclose that the flow leaving such 
(Theoretically, this 
reasoning can be extended to an indefinite number of additional 
shock surfaces. ) 

Unfortunately, the presence of more than one shock surface in a 
flow field introduces complications that can only be handled 
approximately by the present theory. Intersections of shock 
waves produce vortex layers with discontinuities in velocity and 
entropy in the disturbed region behind the shocks. Such dis- 
continuities are impossible to handle by the linearized expres- 
sions, Eqs. (1), because the x and y coordinates of the vortex 
layers vary with ¢ for a ¥ 0. 
then is necessary but is not a serious limitation except where a 


The neglect of any discontinuities 
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secondary shock* is strong. In application, if the vortex layers 
have been found to be of small consequence in the axially sym- 
metric solution, the inclined flow theory may be successfully em 
ployed. 

Numerical-graphical application of the method of character 
istics to a body with multiple shock surfaces is similar to the 
method given in reference 1. Important exceptions, however, are 
the neglect of entropy gradients due to secondary shocks and the 
use of Eq. (7c) for computing w at second-shock fronts. Condi- 
tions at shock intersections are handled by means expounded in 
axially symmetric theories with the neglect of the discontinuities 
associated with vortex layers. 

A more complete discussion on the use of the method of char 
acteristics for inclined flow about a body of revolution, coupled 
with numerical examples on both single and double shock-surface 
configurations, is given in reference 2 

* Here and throughout this paper, the secondary or weak shock is sup 
posed to originate somewhere behind the nose of the body The comments 
are equally valid, however, in cases where the shock from the nose can be 


treated as the weak shock 
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T° THE PAPER, ‘The Solution of Elastic Problems by Electronic 
Analog Computation,”’ Jonathan Winson (JOURNAL OF THE 
AERONAUTICAL SCIENCES, Vol. 17, No. 7, p. 385, July, 1950) has 
undertaken the criticism of another method of machine computa- 
tion analysis suitable for aeroelastic problems. This is the 
method employing what he terms the ‘direct component-to- 
component approach to analog computation 
employed as one of the techniques used on the California In 
stitute of Technology electric analog computer for the solution 
of many aeroelastic problems for the aircraft industry in South- 


This has been 


ern California. 

We should like to take exception to 
made by Mr. Winson, who seems to be 
possibilities of these analog techniques 
Mr. Winson’s experience apparently has only been with the 
REAC type of computation, and he has therefore only found it 
possible to solve problems of a limited degree of complexity 

He states, in Section II of his paper, that “the direct com 
ponent-to-component approach” offers serious difficulties ‘‘in 
that a basically new problem requires a major application of 
effort and ingenuity to invent a suitable analog, with no assur 
ance that such an analog can be devised at all.’’ We wish to 
point out that this statement is untrue. There is an existence 
theorem that states that any conservative linear mechanical 
system has its complete analog consisting only of inductors, 
Further, we have developed 


some of the statements 
unfamiliar with the full 
as applied to this field 


capacitors, and transformers 
general methods for the application of these direct analogies to 
missile and airplane structures and have used them for flutter 
and gust-loading analyses of considerable complexity 
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The CalTech electric analog computer has now been employed 
for comprehensive flutter design analyses in the development of 
several new airplanes. In these studies, elastic structures with a 
large number of degrees of freedom have been represented. 
These include the simulation of wings, fuselages, and stabilizers 
as distributed beams. Control surfaces, engine mountings, and 
landing gears can, of course, be included. Aerodynamic forces 
and torques can be applied to all surfaces, whether fixed or mov- 


able. Currently accepted expressions for the effects of finite 
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aspect ratio, sweepback, downwash, and compressibility cay 
handled. 

The two greatest advantages of our type of analysis are 
the coordinates used in the electrical analogy correspond to 4 
actual displacements and rotations of stations in the struety ; 
and that local modifications of the actual structure result 
equally local modifications of the analogy circuit, permittj 
therefore, rapid consideration of modifications as the design 


analysis proceeds 
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